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KEY  TO  SYMBOLS 


A = Helmholtz  free  energy 

= ith  order  reduced  density  coefficient  of  DCFI  model 
a^j  = coefficient  matrix  of  DCEI  model 
b = constant 

C = direct  correlation  function  integral 
c = molecular  direct  correlation  function 
Cp  = heat  capacity  at  constant  pressure 
d = effective  hard  sphere  diameter 
f = state  dependent  function 

g = pair  distribution  function,  state  dependent  function 
H = enthalpy,  spatial  integrals  of  total  correlation  function 
h = total  correlation  function 
_I  = identity  matrix 
k = Boltzmann's  constant 
k- • = binary  parameter 
M = general  property 
m = integer 

N = number  of  moles,  number  of  components 
n = integer,  refractive  index 
P = pressure 
R = gas  constant 
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r = position  vector 

5 = entropy 

s = position  vector 
T = temperature 
U = internal  energy 

u = velocity  of  sound,  intermolecular  potential 

V = volume 

v = V/N,  molar  volume 

v^  = partial  molar  volume  of  component  i 
X = matrix  of  mole  fraction 
x = mole  fraction,  position 

Y = density  dependent  function  of  hard  sphere  diameter 
Oip  = isobaric  expansivity 

= activity  coefficient  of  component  i 
yv  = thermal  pressure  coefficient 

6 = Kronecker  delta 
£ = parameter  set 

0^  = parameter  of  molecule  i 
0 = e/kT 

e = energy  parameter,  dielectric  constant 

r|  = packing  fraction 

Hy  = isothermal  compressibility 

p = chemical  potential 

v = Stochiometric  coefficient 

p = density 

a = collision  diameter 
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T = inverse  of  reduced  temperature 
cp^  = volume  fraction  of  component  i 
ft  = orientation  normalization  factor 
oo  = angular  orientation  coordinate 

Subscripts 

c = critical  property 

i,j,k  = component 

ij»jk»ik  = pairs  of  components 

m = mixture 

o = reference  property,  saturation  property 
hs  = hard  sphere 


Superscripts 

E = excess  property 
r = reference  state 
* = characteristic  parameter 
o = standard  state 
' = first  composition  derivative 
" = second  composition  derivative 


(as  in  X) 

(as  in  v ) 

~ (as  in  p) 

o (as  in  v ? ) 
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Special  Symbols 
Matrix  quantity 
Partial  molar  property 
reduced  property 
pure  component  property 
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ensemble  average 

angular  orientation  average 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

THERMODYNAMIC  PROPERTIES  OF  COMPRESSED  LIQUIDS  AND 
LIQUID  MIXTURES  FROM  FLUCTUATION  SOLUTION  THEORY 

By 

Yung-Hui  Huang 
August  1986 

Chairman:  John  P.  O'Connell 

Major  Department:  Chemical  Engineering 

For  every  pure  liquid  there  exists  a characteristic 
volume  for  which  the  reduced  bulk  modulus,  v/hjRT= 

(3  (P/RT )/9 p ) j , is  independent  of  temperature.  In  addition 
to  the  possible  physical  significance  of  this  effect,  the 
crossover  point  serves  as  the  basis  for  a correlation  of 
the  spatial  integrals  of  the  statistical  mechanical  direct 
correlation  function  (DCFI)  which  can  be  related  to  the 
density  derivatives  of  thermodynamic  properties  in  the 
fluctuation  solution  theory. 

The  pressure  variation  of  the  volume  and  reduced  bulk 
modulus  of  essentially  all  pure  and  mixed  liquids  are 
correlated  with  a three-parameter  corresponding  states 
correlation  over  temperatures  from  the  triple  to  nearly 
the  critical  point  and  densities  from  the  saturation  to 


xi 


the  freezing  point.  The  parameterization  is  based  on  the 
detailed  analysis  of  the  most  accurate  experimental  data. 

The  three  parameters  each  serve  a particular  mathematical 
purpose,  ensuring  flexible  and  tractable  numerical  calcula- 
tions. 

The  correlation  method  has  been  applied  to  over  200 
pure  and  mixed  liquids  with  errors  probably  within  those 
of  the  experiments,  making  it  significantly  more  general 
and  accurate  than  previous  models.  The  systems  investi- 
gated cover  all  orders  of  complexity  of  molecular  struc- 
ture and  intermolecular  interaction.  These  include 
liquefied  noble  gases,  short-  and  long-chain  hydrocarbons, 
small  and  large  globular  chain  or  ring  compounds,  polar 
and  highly  associating  species,  metallic  liquids,  molten 
salts,  and  polymer  melts.  Mixtures  of  simple  and  complex 
systems  are  carefully  examined. 

The  results  are  relatively  insensitive  to  the  charac- 
teristic temperature  and  all  parameters  may  be  estimated 
by  a group  contribution  method.  Simple  mixing  rules  work 
successfully  in  a one-fluid  form  for  mixtures  with  small 
excess  volume.  In  other  cases,  a single  binary  parameter  is 
adequate  for  good  agreement.  These  findings  indicate  that 
the  direct  correlation  function  integrals  are  less 
sensitive  to  the  exact  details  of  the  intermolecular 
forces  for  dense  liquids  than  other  statistical 
mechanical  quantities. 
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This  correlation,  along  with  accurate  saturation 
values,  will  provide  all  thermodynamic  properties  of 
compressed  liquids  and  liquid  mixtures. 


CHAPTER  1 
INTRODUCTION 


High  pressure  liquids  and  liquid  mixtures  are  pres- 
ent in  many  natural  and  man-made  situations. ^ They  vary 
from  geological  systems  to  high  pressure  processes  and 
techniques.  Thermodynamic  properties  of  compressed 
liquids  and  their  mixtures  are  important  both  in  basic 
scientific  research  and  practical  engineering  applica- 
tions. Thus  the  theoretical  study  of  the  volumetric  beha- 
vior of  compressed  liquids  is  motivated  by  the  concept 
that  the  structure  of  dense  liquids  is  dominated  by  the 

short-range  repulsive  forces  and  that  their  physical  pro- 

2 

perties  can  be  treated  in  terms  of  a hard  sphere  model. 
From  the  practical  point  of  view,  liquid  volumes  are  used 
exclusively  in  equipment  sizing  and  descriptions  of  occu- 
pancy in  natural  reservoirs  and  formations. 

Since  the  pioneering  work  of  Bridgman,*^  there  have 
been  many  volumetric  measurements  of  compressed  pure 
liquids,  whereas  much  less  effort  has  been  dedicated  to 
compressed  liquid  mixtures.  A recent  survey^  of  PVT 
properties  of  liquid  and  liquid  mixtures  listed  about  350 
pure  components  and  170  binary  mixtures.  However,  half  of 
the  numbers  cited  are  restricted  to  very  limited  pressure 
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or  temperature  ranges,  including  many  of  atmospheric  or 
single  isotherm  measurements.  Compression  data  of  pure 

liquids  have  been  correlated  with  equations  such  as  those 

5 7 fl 

of  Tait,  Hudleston,  Hayward,  Chaudhuri,  and  others. 

Despite  their  substantial  value  for  smoothing,  interpola- 
tion, and  computation  of  compressional  thermodynamic  func- 
tions, these  equations  share  common  disadvantages.  First, 
they  contain  no  explicit  temperature  dependence,  so  exten- 
sions by  using  temperature-dependent  parameters  are  still 
uncertain.  Second,  two  parameters  are  insufficient  for  wide 
pressure  range  isotherms,  so  extrapolations  to  either  low  or 
extreme  high  pressures  are  inaccurate.  Third,  their  empiri- 
cal nature  renders  molecular  interpretation  difficult. 

Fourth,  these  equations  do  not  lead  to  mixing  rules  for 
compressed  liquid  mixtures. 

An  acceptable  equation  of  state,  apart  from  quanti- 
tative success,  should  provide  some  insight  into  the 
molecular  characteristics  in  terms  of  a minimum  possible 
set  of  parameters.  Agreement  between  theory  and  experi- 
ment not  only  enhances  one's  faith  in  the  molecular  theory 
but  also  leads  to  the  ability  to  generate  the  volumetric 
properties  with  confidence  from  very  limited  experimental 
efforts.  Theoretical  equations  of  state  which  are  based 
on  the  modification  of  Prigogine's  cell  theory^  such  as 
Flory  et  al.,^  Simha  and  Somcynsky,  and  others  have  been 
applied  to  compressed  liquids  and  liquid  mixtures  using 
scaling  parameters  evaluated  from  low-pressure  density  data.- 
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However,  as  noted  by  Bridgman,^  experimental  behavior  at  low 
pressures  is  a poor  guide  to  what  is  likely  to  happen  at 
high  pressure.  As  a result,  these  equations  of  state  can 
not  correctly  represent  the  PVT  behavior  of  compressed 
liquids  at  high  pressure.  Empirical  modifications  by 
Zoller,^  and  Simha  and  his  coworkers^"^  using  elevated 
pressure  data  make  only  marginal  improvement.  Apparently, 
the  theory  cannot  furnish  an  ideal  representation  of 
experimental  PVT  data  at  elevated  pressure. 

Recently,  the  Tait  equation  has  been  extended  and 
generalized  by  Thomson  et  at.^  This  correlation,  along 
with  the  correlation  developed  by  Hankinson  and  Thomson^ 
for  saturated  liquid  densities,  comprises  the  COSTALD 
method.  These  authors  claim  that  COSTALD  is  the  most 
general  and  accurate  compressed  liquid  density  correlation 
yet  published.  This  is  true  when  compared  with  those 
methods  commonly  used  in  chemical  process  calculations 
such  as  Yen-Woods,  Chueh-Prausnitz , Lee-Kesler,  and 

others.  However,  the  correlation  method  is  purely  empiri- 
cal; application  is  limited  to  nonpolar  and  slightly  polar 
liquids  and  pressures  no  higher  than  700  bar. 

A judicious  alternative  is  to  combine  the  simplicity 
of  an  empirical  equation  with  the  advantage  of  statistical 
mechanical  theory  in  developing  accurate  and  general  rela- 
tions for  thermodynamic  computations.  Fluctuation 
solution  theory  is  explored  here  for  this  purpose.  The 
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basis  of  the  method  is  from  the  statistical  mechanical 
grand  canonical  ensemble  which  relates  composition 
derivatives  of  the  pressure  and  the  chemical  potential, 
thus  the  activity  coefficient,  to  integrals  of  the 
molecular  direct  correlation  function.  As  recognized  by 
O'Connell , 22-23  ancj  Gubbins  and  O'Connell,  ^ the  density 
dependence  of  the  reduced  bulk  modulus  behaves  differently 
and  quite  simply  in  liquid  phase  compared  to  the  gas  phase, 
regardless  of  the  complexity  of  the  intermolecular  forces. 
Hence,  this  work  focuses  on  describing  the  liquid  state 
alone . 

The  most  surprising  feature  of  the  reduced  bulk  modu- 
lus versus  reduced  density  plot  is  that  there  exists  a 
unique  point  of  crossover  of  the  isotherms  for  all  of  the 
compressed  liquids  studied.  This  finding  is  obviously  a 
most  important  contribution  to  any  liquid  state  theory.  In 
view  of  the  well-known  fact  that  it  is  density,  not  pres- 
sure, that  is  the  key  variable  for  compressed  liquids,  a 
corresponding  states  correlation,  based  on  the  density  and 
temperature  expansion  of  reduced  bulk  modulus,  has  been 
constructed  to  predict  accurately  the  pressure  effect  on  the 
compressional  thermodynamic  properties. 

In  the  chapters  that  follow,  detailed  descriptions 
of  the  theoretical  background,  experimental  observations, 
model  development,  and  successful  applications,  including 
pure  and  mixed  liquids  are  presented.  Chapter  2 
describes  the  fluctuation  solution  theory  and  how  useful 
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thermodynamic  properties  are  related  to  the  statistical 
mechanical  direct  correlation  function  integrals. 

Chapter  3 displays  the  reduced  bulk  modulus  for  many  com- 
pressed liquids,  providing  significant  features  that  are 
further  exploited.  For  compressed  liquid  mixtures, 
volumetric  data  are  fully  analyzed  and  deeper  insights  con- 
cerned with  the  excess  properties  are  carefully  evaluated. 
Chapter  4 is  devoted  to  the  process  of  model  construction. 

It  is  demonstrated  how  the  theoretical  concept  can  be  trans- 
formed into  an  actual  working  equation  that  simulates,  to 
within  experimental  accuracy,  the  compressed  liquid  beha- 
vior. Chapter  5 applies  the  model  to  pure  liquids.  Also 
basic  thermodynamic  properties  such  as  entropy,  internal 
energy,  and  others  are  derived  through  the  Maxwell  rela- 
tions. Characteristic  parameters  and  comparison  with  a 
variety  of  pure  liquids  compression  data  are  tabulated.  An 
initial  attempt  at  a group  contribution  correlation  for  the 
parameters  is  also  given.  Chapter  6 extends  the  model  to 
liquid  mixtures.  Mixture  properties  such  as  partial  molar 
volume,  activity  coefficient,  and  others  are  expressed  in  a 
one-fluid  form.  Binary  parameters  and  comparison  with 
liquid  mixture  compression  data  are  tabulated.  Chapter  7 
closes  this  work  with  some  remarks  on  the  numerical  techni- 
ques and  our  present  concept  of  compressed  liquid 
properties . 


CHAPTER  2 

FLUCTUATION  SOLUTION  THEORY  FOR  COMPRESSED  LIQUIDS 

AND  LIQUID  MIXTURES 

Introduction 

Over  the  years  there  have  been  a growing  number  of 
uses  of  statistical  mechanics  in  physical  property 
research.  There  are  two  general  approaches  in  the  devel- 
opment of  liquid  state  theory.  One  may  be  called  the 
formal  approach,  based  on  distribution  functions,  while 
the  other  may  be  called  the  model  approach,  based  on 
partition  functions.  The  former  was  pioneered  by  Mayer 
and  Mayer,  J Kirwood,  ° and  others.  Thermodynamic  proper- 
ties can  be  calculated  through  the  energy  equation  and  the 
pressure  equation  with  known  pairwise  additive  potentials 
and  pair  distribution  functions.  The  latter  approach 
visualizes  a liquid  with  a certain  physical  picture,  such 
as  made  up  of  cells,  holes,  or  molecular  arrangements. 
Thermodynamic  properties  then  can  be  obtained  by  trans- 
forming these  physical  models  into  a partition  function. 
All  suffer  the  disadvantage  of  assuming  pairwise  additi- 
vity of  i nt er m ol ecul ar  potentials. 

Another  alternative  may  be  called  the  fluctuation 
solution  theory. ^ The  compressibility  equation  which 
relates  concentration  derivatives  of  pressure  and 
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chemical  potential  to  the  spatial  integrals  of  the  total 
correlation  function,  and  hence  to  direct  correlation 
function  integrals  (DCFI)  in  the  grand  canonical  ensemble 
is  free  from  this  assumption.  Furthermore,  the  DCFI 
appear  to  be  insensitive  to  the  exact  details  of  inter- 
molecular  forces.  Therefore,  from  a numerical  point  of 
view,  the  theory  involves  integration  of  a simple  model, 
instead  of  taking  the  derivative  of  a complicated  expres- 
sion as  in  the  others.  Experience  shows  that  this  should 
be  less  sensitive  to  errors  of  modeling.  It  is  the  number 
density  and  temperature,  not  the  usual  mole  fraction, 
pressure,  and  temperatures  that  are  the  most  appropriate 
set  of  independent  variables  for  theoretical  work.  How- 
ever, the  equations  allow  for  consistent  relationships  to 
be  made  among  the  variables. 

Fluctuation  Solution  Theory  and 
Direct  Correlation  Function  Integrals 

In  the  statistical  mechanical  grand  canonical 
ensemble,  a relation  exists  between  composition  fluctua- 
tions and  derivatives  of  the  mean  composition  with  respect 
to  chemical  potentials. 


9<Nj> 

'SUi/kT 


T,  V,  pk;I<i 


- <NiXN  > 


<NiNj> 


(2-1) 
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where  the  brackets  denote  a grand  canonical  ensemble 
average.  These  averages  are  also  related  to  the  spatial 
integrals  of  the  radial  distribution  function  by 


<NiXN,>  f 

<NiNj>  - 6ij<Ni>  = J gij(r)dp 

(2-2) 

where  is  the  Kronecker  delta.  As  recognized  by 

Kirkwood  and  Buff,-^  g^j  (r)  is  not  necessarily  limited 
to  be  the  radial  distribution  functions  for  spherical 
molecules;  it  can  also  be  interpreted  as  a more  general 
pair  distribution  function  with  all  possible  orientations. 
Using  the  definition  of  total  correlation  function, 


(2-3) 


A combination  of  equations  (2-1),  (2-2),  and  (2-3)  yields 


9 <N  i > 

J 

Sbi/kT 


T>  V,pk?<i 


5ij<Ni> 


<Ni><Nj> 


<hij(r»u> 


dr 


V 


(2-4) 
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where 


<hi j (r  )> 


0)  = 


h ± j ( rw^co  j )dcoidwj 


= the  angle  averaged  total  correlation  function. 


ft  = J dco  = the  nomalization  constant  of  the 
orientation  dependence. 


r = - r j = the  translational  change  of  the 

position  vector  i and  j. 


Equation  (2-4)  can  be  written  as 


3Pj 


SUj/RT 


T’ V,yk*i 


= 5i jPi 


PiPj 


<hij(r)>0) 


dr 


(2-5) 


This  is  the  basic  equation  of  the  Kirkwood-Buff  solution 
theory.  From  Orstein  and  Zernike,^-*-  the  direct  cor- 
relation function  is  defined  as 


Cij(?i?ja3ia)j)  = hij(ripj0)ia)j)  - Z pk 

JJ  Cik  (-^irka3ia)i<  ) hk  j ^ rkr  j ^ drkdajk 

(2-6) 


The  second  term  on  the  right-hand  side  represents  the 


indirect  correlations  due  to  all  possible  third 
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molecules.  Similarly,  the  angle  averaged  version  of 
Equation  (2-6)  is 


<Cij(r)>w  = <h  — (r)>w  - Z pk 


<Cik(s)>w  <hkj  (r“s)>a)  ds 


(2-7) 


where  <Cij(s)>u  = J J cL  j ( ru^to  j ) da^dio  j 


the  angle  averaged  direct  correlation 
function . 


s = r^  - rk  = the  translational  change  of  position 
vector  i and  k. 

r-s  = r^  - Tj  = the  translational  change  of  position 
vector  k and  j . 

The  angle  averaging  operation  in  the  last  term  of  Equation 
(2-7)  may  not  be  appropriate  for  long-range  potentials.^ 
Integration  of  Equation  (2-7)  yields 


<Cij(r»W 


dr  = p 


<hij(r)>u) 


dr  - p‘ 


Z 

k 


<cik(S)> 


0) 


<hkj(r-s)>w  dr  ds 


(2-8) 
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In  matrix  notation,  Equation  (2-8)  can  be  expressed  as 

C = H - C X H ' (2-9) 

where  X is  a diagonal  matrix  whose  nonzero  elements  are 
mole  fractions.  After  basic  matrix  manipulations,  Equa- 
tion (2-9)  can  be  written  as 

I + H X = (I  - C X)"1  (2-10) 


or 


6ij  + Pi  J <hij^)o)>d?  = <5ij  - Pi 


<cij(?)a,>d?)-1 

(2-11) 


Comparing  Equation  (2-11)  with  Equation  (2-5)  and 
recognizing  that  the  elements  of  the  left-hand  side  of 
Equation  (2-5)  are  the  matrix  inverse  of  a more  desirable 
partial  derivative, 


SPi/RT 


T,pk*j 


(2-12) 


where  C^j  = p j <C^j(r)>w  dr  = the  direct  correlation 
function  integrals.  This  is  the  basic  connection  between 
thermodynamic  variables  and  direct  correlation  function 
integrals. 
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Thermodynamic  Consequences 

From  the  isothermal  Gibbs-Duhem  equation,  the  compo- 
sition derivative  of  the  pressure  and  the  chemical  poten- 
tials are  related  to  each  other  by 


9P/RT 


T,pk*j 


= Z 
i 


3yi/RT 

Tp~ 


T,pk*j 


(2-13) 


Substitution  of  Equation  (2-12)  into  Equation  (2-13) 
yields 


9P/RT 

Tp- 


T,P 


k*j 


= Z 
i 


:i  ( 1-Ci  j ) 


( 2-14a ) 


FT 


(2-14b) 


where  the  subscript  m indicates  mixture  properties  and 
x-p  is  the  isothermal  compressibility.  Summing  all  equa- 
tions (2-14)  yields 


9P/RT 

“Tp— 


T,  x 


= Z Z X^d-Cy.-) 
i J J J 


(2-15a) 


1 

pmwTmRT 


( 2-15b ) 
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The  partial  molar  volume  V j then  can  be  expressed  in 
terms  of  density  and  the  DCFI 


V . 
J 


Z x 


i(l  c j.  j ) 


m 


l Z xix-(l-CijJ 


i J 


(2-16) 


where  pm  is  mixture  molar  density.  An  isothermal  change 
of  pressure  can  be  obtained  by  integration  of  Equation 
(2-14) 


P-Pr 

RT 


(p-pr)  - Z 


p.C.  . 
i ij 


T,P 


dp 


i*j 


(2-17) 


while  an  isothermal  change  of  chemical  potential  is 


T,Pk*j 


(2-18) 


where  the  superscript  r denotes  a reference  property. 

This  can  be  any  fluid  state,  including  ideal  gas, 
saturation,  or  other  well-determined  conditions. 
Naturally,  the  model  for  the  must  be  appropriate  over 

the  entire  range  of  the  integration.  In  phase  equili- 
brium calculations,  it  is  usually  more  convenient  to  work 
with  activity  coefficient  than  chemical  potential. 
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Ui(T»p)  = U°(T,pr)  + RTlnx^iCT,  p) 


(2-19) 


where  the  superscript  o denotes  the  standard  state  pro- 
perties. Substitution  of  Equation  (2-19)  into  Equation 
(2-12)  yields 


T,pk*j 


1-C 


i j 


P 


(2-20) 


Composition  integration  of  equation  (2-20)  yields 


In 


Z 

j 


1-C.  . 

(— rM)  dp  ■ 


For  a pure  component, 


(2-21) 


3P/RT 

“5p~ 


T 


1-C 


(2-22a) 


1 ( 2-22b ) 

PHjRT 

This  is  the  reduced  bulk  modulus  of  pure  liquids,  the 
quantity  of  primary  interest  in  the  present  work. 


CHAPTER  3 

VOLUMETRIC  BEHAVIOR  OF  COMPRESSED  LIQUIDS 
AND  LIQUID  MIXTURES 


Introduction 

Theoretical  investigations  of  the  physical  properties 
of  dense  liquids  indicate  that  the  structure  of  liquids, 
static  or  dynamic,  is  similar  to  that  of  a rigid-body 
system.  In  order  to  examine  the  effects  of  the  differ- 
ences between  real  and  rigid-body  systems,  the  temperature 
and  density  effects  should  be  separated  requiring  experi- 
mental studies  at  high  pressures.  In  describing  the  ther- 
modynamic state  of  a liquid,  the  applied  pressure  is 
considered  high  if  it  is  comparable  to  or  greater  than  the 
kinetic  pressure.^  The  thermodynamic  equation  of  state 
is 


p - t ( ^ ^ ^ ) 

P " 1 ^TTV  [JVJT 


(3-1) 


3 P 

where  the  thermal  pressure  coefficient , T(-^y)y , is  the 
kinetic  pressure  resulting  from  the  molecular  motion  and 
the  energy-volume  coef f icient , (yy) y , is  a measure  of  the 
internal  pressure  arising  from  the  intermolecular  forces. 
When  the  energy-volume  coefficient  (yy)j  is  zero,  the 
external  pressure  P is  equal  to  the  kinetic  pressure 

T(It)V  ’ and 
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(|t)v  = T (3-2^ 

This  is  the  case  for  a rigid-body  fluid.  When  the  applied 
pressure  P is  close  to  the  kinetic  pressure,  a perturba- 
tion approach  is  suggested  using  a rigid  body  reference 
with  a small  perturbation  attractive  potential  term. 

Under  these  circumstances,  the  volume  plays  the  dominant 
role  in  affecting  the  molecular  behavior  in  liquids, 
because  the  molecules  are  so  closely  packed  together  that 
the  forces  between  them  arise  almost  exclusively  from  the 
short-range  repulsive  forces. 

There  has  been  much  theoretical  speculation  concern- 
ing this  physical  picture  of  liquids.  The  first  contro- 
versy comes  from  the  experimental  evidence  that  the 
internal  energy  passes  through  a minimum  which  is  charac- 
teristic of  the  liquid  and  a function  of  temperature.3^ 

This  brings  in  a temperature  dependent  effective  hard-core 
diameter  which  would  be  necessary  to  make  Equation  (3-2) 
behave  like  Equation  (3-1).  Further  experimental 
evidences  for  highly  compressed  liquids  indicate  that 
theoretical  values  of  the  hard  sphere  liquid  deviate 
systematically  from  the  experimental  values.^  This 
suggests  that  the  effective  hard-core  diameter  becomes 
density  dependent  at  high  packing  fractions  (pa^>0.93)  due 
to  the  softness  of  the  repulsive  potential  and/or  orienta- 
tional ordering  in  compressed  liquids.  Figure  3-1  shows  the 
reduced  bulk  modulus  of  a hard  sphere  liquid  calculated 
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Figure  3-1.  Theoretical  DCFI  of  Effective  Hard  Sphere 
Liquid 
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from  the  Percus-Yevick  compressibility  equation.  The 
Verlet-Weis  algorithm  is  adopted  to  simulate  the  tempera- 
ture and  density  dependent  effective  hard-core  diameters 
of  Weeks-Chandler-Andersen  perturbation  theory.2  Detailed 
calculations  are  summarized  in  Appendix  A.  Figure  3-2 
shows  the  experimental  reduced  bulk  modulus  of  liquid 
argon  calculated  from  Twu's  algorithm.-56  Comparison  of 
Figures  3-1  and  3-2  shows  that  the  differences  are  signi- 
ficant, especially  at  high  pressure  where  the  crossover  is 
observed  for  argon  and  the  temperature  dependence  is  very 
small . 

Since  theoretical  calculations  of  compressed  liquids 
are  improbable,  the  alternative  is  to  analyze  fully  the 
experimental  data  available  and  develop  an  empirical  rela- 
tionship. There  are  several  interesting  findings  concern- 
ing the  volumetric  behavior  of  compressed  liquids  and 
liquid  mixtures.  The  earliest  one  is  made  by  Bridgman  on 
the  pressure  effect  of  isobaric  expansivity.-5  At  low 
pressure,  the  thermal  expansion  coefficient  of  liquids 
increases  with  rising  temperature.  At  high  pressure,  this 
relation  is  reversed.  The  reversal  in  the  sign  of  the 
temperature  derivative  of  isobaric  expansivity  is  due  to 
the  nonlinear  effect  of  the  intermolecular  forces.  It  is 
interesting  to  relate  this  phenomenon  to  the  behavior  of 
constant  pressure  heat  capacity.-5 
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pvc 


Figure  3-2.  Experimental  DCFI  of  Liquid  Argon 
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3CP, 

yp-;T 


( 3-3a ) 


= -VT  (a 


2 daP, 

p +ir) 


( 3 -3b ) 


where  ap  is  the  isobaric  expansivity.  The  effect  found 
by  Bridgman  raised  the  possibility  of  a minimum  in  Cp  as  a 
function  of  pressure. 

Similar  observations  on  the  change  of  internal  energy 
with  pressure  were  also  made  by  Bridgman.^ 


_ y(3V\  (3-4) 

T ’ UTTJP  ^TF't  u 

The  pressure  at  which  the  left-hand  side  of  (3-4)  is  zero 
is  given  by  the  same  form  as  Equation  (3-2). 


P 


-T  ( 3 ^ ) 

U7T;P 

-TV— 

^TFJT 


= T ( 


3 P 


FT'  V 


), 


(3-5) 


These  variations  were  further  explored  by  Jenner  and 
his  coworker^"^  on  the  compression  measurement  of 
1-bromoalkanes  from  C£  to  C -j  in  the  liquid  state.  The 
applied  pressure  is  up  to  6 Kbar  in  the  temperature  range 
between  203°K  and  448°K.  For  all  of  the  bromides 
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investigated,  there  exists  a well-determined  pressure  for 
which  the  isobaric  expansivity  is  independent  of  tempera- 
ture. This  inversion  pressure  is  further  confirmed  by 
experimental  data  from  the  piezothermal  method  of  Ter 


dioxide,  and  n-butane.  The  curves  of  the  pressure  deriva- 
tives of  internal  energy  and  of  heat  capacity  show  similar 
behavior  as  a function  of  pressure.  In  addition  to  the 
observed  minimum  at  different  temperatures,  all  of  the 
isotherms  pass  through  a single  point  at  the  inversion 
pressure  noted  above.  Further,  the  isothermal  variation 
of  the  heat  capacity  difference  Cp-Cy  exhibits  the  same 
crossover  at  the  inversion  pressure. 

In  all  of  these  cases,  there  must  exist  a thermodyna- 
mic relationship  at  the  inversion  pressure  p^.  First, 


Minassian  and  his  coworkers^1’42  on  liquid  water,  carbon 


daPi  1/32V 


2 


0 


(3-6) 


leading  to 


Pi 


(3-7) 


Second , 
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-T  ( 


_afV) 

3T2  Pi 


(3-8) 


where  by  is  a constant.  Combination  of  Equation  (3-7)  and 
(3-8)  yields 


TV 


(3-9) 


Third , 


C D ~ 


C - T ( ^ P 'i 

CV  - T(TT)Pi(FT)V 


= JVaPx  t\i 

c. 


apiYV 


(3-10) 


where  b2  is  another  constant.  From  Equation  (3-10),  the 
thermal  pressure  coefficient  yy  can  be  written  as 


Y 


* 

V 


FJ  “Pi 


(3-11) 


From  the  triple  product  rule,  the  isothermal  compressibi- 
lity Hy  can  then  be  given  as 
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h?  = 


aPi 


*1 


(3-12) 


This  is  true  only  when  inversion  pressures  are  exactly  the 
same  for  all  of  the  properties  investigated  above.  Unfor- 
tunately, data  analysis  on  the  isothermal  compressibility 
is  unavailable  in  Jenner's  work  to  confirm  this  com- 
pletely. Careful  analysis  of  this  property  will  be  des- 
cribed in  next  section. 


Pure  Liquids 

Volumetric  properties  of  large  globular  and  long-chain 
hydrocarbons  have  been  measured  by  Grindley^  for  tempera- 
tures from  298°K  to  453°K  and  pressures  from  zero  to 
8 Kbar.  The  data  are  adequate  to  determine  the  first  and 
second  volume  derivatives  of  the  internal  energy  and 
entropy.  The  first  isothermal  volume  derivatives  yield 
the  thermodynamic  equation  of  state. 

dt/3Ss  /8U> 

p = - Wi  <3-i3) 

Both  terms  on  the  right-hand  side  of  (3-13)  are  tempera- 
ture dependent.  However,  scaling  can  be  achieved  by  a 
two-parameter  corresponding  states  relation.  The  second 
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isothermal  volume  derivatives  are  related  to  the  bulk 
modulus 


In  contrast  to  the  external  pressure  expression  (3-13), 
where  the  entropy  term  dominates  throughout  the  whole 
density  region,  the  bulk  modulus  has  a large  contribution 
from  both  terms.  The  entropy  term  is  larger  at  low  densi- 
ties while  the  energy  term  is  larger  at  high  densities. 
This  can  be  seen  clearly  from  Figure  3-3  for  the  second 
isothermal  volume  derivative  properties  of  n-nonane.  The 
crossing  of  the  bulk  modulus  isotherms  is  due  to  the 
additive  effect  of  both  entropy  and  energy  term,  though 
crossing  of  the  energy  term  also  appears  at  higher 
densities . 

The  thermodynamic  implications  of  this  behavior  have 
been  analyzed  further  in  this  work.  Our  examination  of 
the  experimental  reduced  bulk  modulus  data  available 
reveals  a similar  feature:  A density  at  which  the  reduced 

bulk  modulus,  which  is  related  to  the  integral  of  the 
direct  correlation  function  of  Chapter  2,  is  independent 
of  temperature.  Figure  3-4  shows  the  experimental  beha- 
vior of  the  DCFI  for  methane  from  Mollerup's  computer 
programs.^  At  the  critical  point,  the  DCFI  is  unity 
because  the  reduced  bulk  modulus  is  zero.  At  the  ideal 
gas  limit,  the  DCFI  is  zero,  and  hence  the  reduced  bulk 


(3-14) 


MOLES /cm 
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Figure  3-3.  Entropy  and  Energy  Contributions  to 
the  Bulk  Modulus  of  n-nonane 
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PVC 

Figure  3-4.  Experimental  DCFI  of  Liquid  Methane 
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modulus  is  unity.  Generally,  the  DCFI  is  small  and  posi- 
tive (less  than  unity)  in  the  vapor  phase;  it  is  large  and 
negative  in  the  liquid  phase.  The  high  density  limits  of 
the  isotherms  are  the  freezing  line,  while  the  low  density 
limits  of  the  liquid  isotherms  are  the  saturated  liquid. 
Figures  3-5  and  3-6  show  the  DCFI  of  two  sets  of  isomers: 
n-nonane,  3 , 3-diethylpentane,  and  n-heptadecane,  5,5- 
dibutylnonane  in  the  liquid  state  alone.  All  the  experi- 
mental curves  behave  similarly  with  relative  insensiti- 
vity to  the  temperature  in  the  higher  density  range,  and 
the  crossover  is  apparent.  By  shifting  the  curves  to 
superimpose  the  crossover  DCFI  and  density,  complete  over- 
lap can  be  obtained  by  rotating  the  graphs  slightly  to 
bring  the  isotherms  into  coincidence.  This  provides  the 
basis  for  a three-parameter  corresponding  state  correla- 
tion of  the  reduced  bulk  modulus.  Similar  behavior  has 
been  found  for  essentially  every  other  substance  studied 
except  water.  Figure  3-7  shows  the  DCFI  of  water  generated 
from  the  computer  programs  of  the  National  Bureau  of 
Standards. ^ The  temperatures  are  from  boiling  to  the 
critical  point  and  pressures  up  to  10  Kbars.  Except  for 
the  low  temperature  isotherms  (below  50°C),  water  exhibits 
the  same  volumetric  behavior  as  the  other  substances  over 
a very  broad  range  of  state  conditions.  The  differences 
at  low  temperature  can  be  attributed  to  hydrogen  bonding 
which  decreases  rapidly  with  temperature,  a conclusion 
which  has  been  reached  on  evidence  from  diverse  sources.^ 
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Figure  3-5. 


pvc 

Experimental  DCFI  of  Liquid  n-nonane  and 
n-heptadecane 
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Figure  3-6.  Experimental  DCFI  of  3 ,3-diethylpentane  and 
5,5-dibutylnonane 
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Volume  and  temperature  have  been  chosen  as  the 
independent  variables  since  they  are  the  most  suitable  for 
theoretical  treatment.  In  general,  volume  scaling  is  the 
result  of  an  increase  in  molecular  size,  while  temperature 
scaling  is  due  to  the  ratio  of  potential  to  kinetic  energy 
effects  and  variation  in  rotational  degrees  of  freedom. 

The  relative  insensitivity  of  the  DCFI  with  temperature  in 
the  compressed  state  indicates  that  internal  degrees  of 
freedom  are  not  appreciably  altered  when  the  density 
changes  above  twice  the  critical  density.  The  crossover 
reduced  bulk  modulus  variation  which  leads  to  DCFI  scaling 
arises  probably  from  the  combined  effects  of  size  and 
inter  molecular  forces  because  of  the  additive  effects  of 
entropy  and  energy  contributions  to  the  crossing  of  the 
liquid  bulk  modulus.  To  better  understand  this  conjec- 
ture, two  separate  studies  from  Gibson  and  Loeffler47  and 
Grindley  and  Lind4®  can  be  considered.  For  the  same  sized 
molecules,  if  a polar  group  is  introduced  as  in  the  former 
case,  or  upon  charging  a set  of  neutral  molecules  to  form 
a molten  salt  as  in  the  latter  case,  the  major  contribu- 
tion to  the  decreased  pressure  and  isothermal  compressi- 


bility (or  increased  DCFI)  comes  from  an  increased  (£iL) 

9 V T 

g 2 u 

which  makes  ( — „ — more  positive.  Only  a small  contribu- 
3 \l 

2 

tion  comes  from  a decrease  of  (|4)_  and  (-- -|) 
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Liquid  Mixtures 


Insights  into  the  thermodynamic  properties  of  liquid 
mixtures  can  be  obtained  from  the  excess  volume  as  a 
function  of  pressure,  temperature  and  composition.  An 


where  Vm  is  the  mixture  volume  and  V°  is  the  volume  of  pure 
component  i.  The  absolute  value  of  excess  volume  is  small 
in  comparison  to  the  mixture  volume  itself,  whereas  the 
variations  of  the  excess  volume  with  pressure,  temperature 
and  composition  can  be  important  and  comparable  to  those  of 
the  mixture  volume. 


equation  of  state  based  on  the  excess  volume  can  be 
written  as 


(3-15) 


TP~  T , x_  Vmx*r 


E 


(3-16) 


m P ,x  = 


(3-17) 


and  for  binary  mixtures, 


( 3-18a ) 


( 3 — 1 8 b ) 
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( 3- 19a ) 


( 3- 1 9b  ) 


Since  very  few  measurements  have  been  made  to  deter- 
mine the  equation  of  state  for  liquid  mixtures,  signifi- 
cant progress  in  their  high  pressure  thermodynamics  has 
not  occurred.  However,  it  is  known  that  the  experimental 
behavior  of  the  excess  volume  as  a function  of  temperature 
and  composition  is  exceedingly  complicated  due  to  a 
variety  of  physical  effects  and  intermolecular  interac- 
tions. One  behavior  of  interest  is  the  pressure  depen- 
dence of  the  excess  volume  of  liquid  mixtures.  Figure  3-8 
shows  the  excess  volumes  of  the  carbon  monoxide-methane 
system  from  Calado  et  al.^  and  the  n-heptane-ethanol 
system  from  Ozawa  et  al.^  The  absolute  value  of 


high  pressure  mixtures  for  a wide  variety  of  liquids 
show  more  ideal  volumes  of  mixing  than  do  low  pressure 
systems.  Generally,  the  pressure  influence  is  very  signi- 
ficant. For  simple  liquid  mixtures,  the  excess  volume  is 
small  and  negative  at  low  pressure,  becoming  almost  zero 
at  high  pressures.  For  n-alkane  mixtures,  the  excess 
volume  is  negative  at  low  pressure,  but  it  eventually 


g y u 

(-^p)  t x rapidly  decreases  with  pressure,  implying  that 
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n-heptane-ethanol  (Ozawa,  et  al.) 


Pressure,  MPA 

CO/CH^  at  xco=0.5  (Calado,  et  al.) 


Figure  3-8.  State  Dependence  of  Excess  Volume 
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becomes  small  and  positive  at  high  pressures.  For 
n-alkane  and  1-alkanol  mixtures,  the  excess  volume  is 
positive  at  low  pressures  becoming  negative  at  high  pres- 
sures. For  aqueous  systems,  the  excess  volume  is  large 
and  negative  at  low  pressures,  becoming  small  and  positive 
at  high  pressures.  These  are  all  closely  related  to  the 
packing  effect  at  highly  compressed  states. 


CHAPTER  4 

CORRESPONDING  STATES  CORRELATION  FOR  THE 
DIRECT  CORRELATION  FUNCTION  INTEGRALS 


Introduction 

The  practical  value  of  the  rigorous  formulas  derived 
in  Chapter  2 and  experimental  observations  discussed  in 
Chapter  3 depends  on  a mathematical  function  to  simulate 
the  DCF  I behavior.  The  ideal  model  for  DCFI  would  yield  a 
robust  mathematical  form,  giving  all  of  the  variations 
nature  shows  with  a minimum  of  experimentally  accessible 
parameters  that  depend  only  weakly  on  thermodynamic 
states.  Treatments  using  the  corresponding  states  princi- 
ple (CSP)  have  this  desirable  feature.  However,  most 
corresponding  states  correlations  are  based  on  the  criti- 
cal properties  and  other  macroscopic  fluid  measurements 
such  as  saturation  properties,  so  they  are  limited  to 
classes  of  substances  that  have  simple  intermolecular 
forces.  Furthermore,  these  properties  are  not  always 
measurable  as  in  the  case  of  heavy  hydrocarbons  and  coal 
liquids.  However,  a corresponding  states  correlation  for 
the  DCFI  would  be  very  useful  in  chemical  process  calcula- 
tions if  suitable  parameters  are  found.  In  this  chapter, 
a three-parameter  corresponding  states  correlation  has 
been  developed  to  describe  the  DCFI  of  all  liquids.  The 
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three  characteristic  parameters  of  the  correlation  are  the 
volume,  V*,  and  DCFI,  C*,  at  the  crossover  of  the  iso- 
therms along  with  a temperature,  T*,  which  is  in  the  range 
of  the  critical  temperature. 

Model  Development 

The  first  apparent  utility  of  the  DCFI  to  describe 
the  relationship  of  pressure  and  density  for  liquids  was 
described  by  Brelvi51  and  Brelvi  and  O'Connell.52-53  They 
noted  that  density,  not  temperature,  was  the  dominant 
variable  and  the  DCFI  of  many  substances  could  be  corre- 
lated in  a one-parameter  corresponding  states  form.  In 
principle,  their  characteristic  volume  parameter  could  be 
determined  from  a single  compression  measurement.  This 
correlation  holds  within  5-10%  for  a variety  of  substances 
over  the  dense  liquid  region.  However,  the  correlation  is 
uncertain  at  lower  densities  where  the  experimental  data 
are  sparse  and  temperature  effect  becomes  significant. 
Further  examination  by  Mathias  and  O'Connell54-55  showed 
that  the  temperature  dependence  could  be  taken  into 
account  with  a characteristic  temperature  in  a two- 
parameter  corresponding  states  form.  Since  temperature 
dependence  was  included,  the  results  at  low  densities  were 
improved.  Unfortunately,  the  correlation  used  hard  sphere 
expres-sions  corrected  by  a linear  density  term  so  appli- 
cations to  high  pressure  were  not  of  high  accuracy. 
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In  the  present  work,  liquid  compression  data  have 
been  analyzed  even  more  carefully  and  the  most  general 
form  of  the  corresponding  states  relation  is  shown  in 
Figure  4-1.  The  behavior  of  the  liquid  DCFI  for  methane 
is  plotted  in  reduced  coordinates  along  with  some  exten- 
sions using  data  for  n-heptadecane.  At  the  crossover 
point,  the  reduced  density  and  DCFI  are  defined  to  be 
unity.  While  variations  away  from  this  characteristic 
density  and  DCFI  are  somewhat  dependent  on  the  tempera- 
ture, they  can  be  correlated  in  reduced  form. 

Model  Parameterization 

The  present  correlation  takes  advantage  of  the 

similarity  of  the  crossover  behavior  that  all  dense  liquid 

DCFI  data  show.  The  success  of  the  model  relies  on  the 

construction  of  an  equation  of  state  that  can  describe  the 

reference  substance  accurately.  Since  this  work  has  been 

concerned  with  formulating  the  thermodynamic  properties  of 

petrochemicals  and  coal-derived  liquids,  methane  is  used 

as  the  reference  substance.  The  basic  form  of  the  cor- 
relation is  chosen  to  express  the  reduced  DCFI,  £=C/C*, 

as  the  polynomial  expansions  of  the  reduced  density, 

p=pV*,  and  the  inverse  reduced  temperature,  x=T*/T. 

m n 

C = Z Z aA  , (T)J_1(p)i_1 
i=l  j=l  J 


(4-1) 
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p/p* 


Figure  4-1. 


Corresponding  States  Correlation  of 
Liquid  DCFI 
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where  a^j  is  the  coefficient  matrix  of  the  reduced 
variables . 

After  a thorough  test  using  the  liquid  DCFI  data  of 
methane  generated  from  Mollerup's  computer  program,  the 
optimal  expression  has  been  determined  as  the  lowest  order 
polynomial  that  can  achieve  the  same  accuracy  as  his  equa- 
tion. In  this  case,  m=4  and  n=3.  The  coefficient  matrix 
a^j  of  equation  (4-1)  is  given  in  Table  4-1.  The  agree- 
ment with  experimental  data  overall  is  better  than  0.2% 
except  for  low  density  points  near  the  critical 
isotherms.  The  valid  range  of  the  reduced  variables  is 
0.7£pV*S1.3  and  0 . 5£  T / T 0 . 9 9 . This  range  could  be 
expanded  by  use  of  a modified  form  of  the  polynomial  and 
data  for  substances  at  lower  reduced  temperatures  and 
higher  reduced  densities.  These  coefficients  give  the 

crossover  value  of  C = 1 to  within  1%  over  the  specified 
temperature  range.  The  reduced  bulk  modulus  then  can 
be  written  as 

9P/RT  _ 1 

5 p T ~ 

= 1 - C*£  (4-2) 


Results  and  Discussions 

Equation  (4-1)  with  the  coefficient  matrix  of  Table 
4-1  has  been  applied  to  all  of  the  substances  for  which 
liquid  DCFI  data  have  been  found.  The  compressibilities 
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TABLE  4-1 


Coefficient  Matrix  a—  for  Equation  (4-1) 


\ 1 ■ 

m\ 

1 

2 

3 

1 

9.8642 

-10.191 

-1.5356 

2 

-28.465 

30.864 

6.0294 

3 

27.542 

-32.898 

-8.7130 

4 

-8.2606 

12.737 

4.0170 

42 


are  generally  from  fitting  volume  changes  with  pressure 
and  sonic  velocities.  The  present  reduced  bulk  modulus 
data  bank  includes  33  substances  ranging  from  noble  gases 
to  petrochemicals  to  water  and  methanol.  (Compressions  of 
many  more  substances  are  discussed  in  Chapter  5.)  Table 
4-2  lists  all  of  the  characteristic  parameters  and  the 
data  ranges  from  which  the  parameters  are  evaluated.  The 
values  of  V*  are  about  30-40%  of  the  critical  volume  while 
the  values  of  T * are  normally  within  20%  of  the  critical 
temperature.  The  values  of  C*  vary  from  -7  for  very  small 
molecules  to  -167  for  large  polyatomics.  The  results  are 
in  excellent  agreement  with  the  experimental  data.  For 
some  substances  with  wide  data  ranges,  the  fitted  pair 
parameters  coincide  with  the  crossover  point  observed 
experimentally.  For  substances  with  lower  pressure 
ranges,  the  fitted  pair  parameters  are  close  to  the  cross- 
over point  obtained  from  graphical  extrapolation  of  the 
experimental  results.  This  is  also  true  when  water  is 
considered  individually  at  high  and  low  temperatures. 

There  are  two  sets  of  characteristic  parameters  for  water 
at  different  temperature  ranges.  They  are  listed  sepa- 
rately because  water  behaves  unusually  at  low  tempera- 
tures56 (below  75 0 C ) . 

Discussion  of  the  results  can  be  developed  according 
to  the  complexity  of  the  substance,  the  range  of  the  state 
conditions,  and  the  accuracy  of  the  experimental  results. 


Characteristic  Parameters  and  Comparisons  with  Liquid  DCFI  Data 
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In  the  case  of  monatomic  substances  from  Streett  et 
ai.56,57  gncj  diatomic  species  from  Bender,-50  the  tempera- 
ture ranges  from  just  above  the  triple  point  to  near  the 
critical  point,  and  pressures  from  saturation  to  freezing 
point.  The  agreement  with  the  present  correlation  is 
within  1%.  The  ratios  of  the  characteristic  parameters  to 
the  critical  properties  are  slightly  lower  than  the  refer- 
ence methane  values,  while  the  characteristic  DCFI  values, 
-C*,  are  slightly  higher.  In  the  case  of  short-chain 
hydrocarbons  from  the  National  Bureau  of  Standards,60-6-5 
the  reduced  temperature  range  is  between  0.5  and  0.97,  and 
the  pressures  are  from  saturation  to  700  bars.  The  agree- 
ment with  the  present  correlation  is  within  1 %.  The  ratio 
of  the  characteristic  volume  to  critical  volume  is  about 
0.35  and  the  ratio  of  the  characteristic  temperature  to 
critical  temperature  is  about  0.92.  In  the  case  of  long- 
chain  hydrocarbons  from  Snyder  and  Winnick,65  the  tem- 
peratures are  25°  to  85°C  and  pressures  from  about  1 bar 
to  over  5 Kbars  or  the  freezing  pressure,  whichever  is 
lower.  The  agreement  with  the  present  correlation  is 
within  2.5%.  This  is  the  worst  for  all  of  the  substances 
tested.  However,  examination  of  the  deviation  patterns 
indicate  that  the  major  source  of  the  disagreement  comes 
from  the  atmospheric  measurements  of  each  isotherm.  There 
are  several  possibilities  related  to  this  observation. 
First,  the  extrapolation  of  high  pressure  data  to  atmos- 
pheric pressure  either  by  a graphic  technique  or  by  a 
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specific  equation  may  cause  errors,  though  it  is  a common 
practice  in  high  pressure  work.  Second,  the  adoption  of 
the  atmospheric  measurements  from  the  other  sources  may 
cause  inconsistency  in  the  complete  data  set.  Third,  the 
initial  compressibility  of  the  long-chain  hydrocarbons  may 
behave  differently  from  that  of  methane.  If  the  devia- 
tions do  result  from  the  experimental  uncertainty  of  the 
atmospheric  DCFI,  then  the  agreement  with  the  present 
correlation  is  definitely  within  the  experimental  errors. 
Similarly,  for  long  chain  and  large  globular  hydrocarbons 
of  Grindly,43  the  temperatures  are  from  25°  to  190°C  and 
pressures  from  zero  to  8 Kbars.  The  agreement  with  the 
present  correlation  is  within  2%.  The  experimental  preci- 
sion of  the  lowest  pressure  compressibilities  is  only  1 or 
2?o,  while  above  400  bars,  the  precision  is  about  .5%. 
Comparisons  of  the  experimental  data  with  the  calculated 
results  indicate  that  the  deviation  patterns  are  consist- 
ent with  the  previous  assertions  that  the  disagreement 
between  the  experimental  DCFI  and  the  present  correlation 
at  zero  pressure  is  due  to  the  experimental  uncertainty  at 
this  state.  For  these  larger  hydrocarbons,  the  ratio  of 
the  characteristic  to  critical  volume  is  about  0.31  and 
the  ratio  of  the  characteristic  to  critical  temperature  is 
between  0.8  and  0.9.  The  only  exception  is  3,3-diethyl- 
pentane  which  has  an  unusually  high  ratio  of  the  charac- 
teristic to  critical  temperature;  this  is  associated  with 
the  abnormally  large  disagreement  of  DCFI  at  zero  pressure 
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where  the  temperature  effect  is  significant.  For  benzene 
and  some  of  its  derivatives,  the  work  of  Gibson  and  Loef- 
is  ranked  as  one  of  the  most  accurate  among  all 
high  pressure  measurements.  The  temperatures  are  from  25° 
to  85 0 C and  pressure  from  zero  to  1 Kbars.  The  agreement 
with  the  present  correlation  is  within  0.5?o  which  is 
significantly  better  than  for  the  rest  of  the  substances. 
This  is  partly  due  to  the  insensitivity  of  the  experimen- 
tal DCF  I with  temperature  at  low  pressures  for  these 
substances.  Figure  4-2  shows  the  experimental  behavior  of 
DCFI  for  benzene.  The  value  of  -C*  increases  from  liquid 
to  liquid  in  the  order  of  benzene,  chlorobenzene,  bromo- 
benzene,  nitrobenzene  and  aniline.  This  is  consistent 
with  the  argument  in  the  previous  chapter  that  the  intro- 
duction of  a polar  group  acts  on  the  liquid  DCFI  of  ben- 
zene in  the  same  way  as  the  increase  of  energy  volume 

coefficient.  It  changes  the  attractive  potential  and 
2 

9 U 

makes  ( more  positive.  The  ratio  of  the  characteris- 

3 \IL  1 

tic  to  critical  volume  is  about  0.34  and  characteristic  to 
critical  temperature  is  about  0.85.  In  the  case  of  the 
molten  salts  from  Grindly,1^^  the  temperatures  are  from  90° 
to  160°C  and  pressures  from  zero  to  5 Kbars.  The  agree- 
ment with  the  present  correlation  is  within  0.8?o.  The 
ratio  of  the  characteristic  to  critical  volume  is  about 
0.32  and  characteristic  to  the  estimated  critical 
temperature  is  about  1.  Comparison  of  -C*  values  for 
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pvc 

Figure  4-2.  Experimental  DCFI  of  Liquid  Benzene 
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tetrabutylammonium  tetrabuty lborate  and  5 , 5-dibuty lnonane 
which  have  nearly  the  same  V*  confirms  that  the  charges 
have  the  same  effect  on  liquid  DCFI  as  a polar  group  does. 
For  highly  associating  species  such  as  water^’^6  and 
methanol^  the  agreement  with  the  present  correlation  over 
a broad  range  of  state  conditions  is  as  good  as  the  other 
substances;  only  lower  temperature  water  is  anomalous. 


CHAPTER  5 

APPLICATION  OF  THE  MODEL  TO  OTHER  PURE  LIQUIDS 


Introduction 

Thermodynamic  properties  of  a pure  liquid  can  be  cal- 
culated from  the  DCFI  model  using  the  parameter  developed 
in  the  last  chapter.  In  addition  to  isothermal  compressi- 
bilities, there  are  compression  data  at  high  pressures,  so 
model  parameters  based  on  the  pressure  equation  and  com- 
pression data  can  be  obtained  for  these  substances  as  well. 
These  parameter  tables  include  almost  all  of  the  signifi- 
cant liquid  compression  data  available.  The  crossover 
parameters  appear  to  be  group  additive  and  related  to 
critical  properties.  Some  of  the  group  parameters  for 
characteristic  volume  V*  and  characteristic  DCFI  C*  are 
tabulated . 


Derived  Thermodynamic  Properties 
Isothermal  integration  of  the  DCFI  model  from  a 
reference  density,  pQ  to  final  density,  p,  yields  the 
pressure,  P 
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P = P0  + P*RT{  (p-p0)-C»[  (9.8642— J.-1.91-.1;?  3 56  ) (p-pQ)- 


(14.2325-1M21-  M^)(p2-^)  + (9.18067-J^9i6-^°^2) 

T T2  T T2 

(p3-p^)-(2.06515-— -1-8425-1.-00^25)  (p4_p4)]} 

T TZ  0 

(5-1) 


where  PQ  is  the  pressure  at  pQ  and  T.  Substitution  of 
the  DCF  I model  into  Equation  (2-22b)  yields  the  liquid  bulk 
modulus  h y ^ , 


k"1  = pRT  {l-C*[ (9.8642- 


10.191  1.5356  ' ,23  /|65  30.864 


T 


^|^)p  + (2  7.54  2-21^j-8-.7^30  )p2_(8.2  60  6-^ilZlZ- 
TZ  T T2  T 


4 . 0170  v ~3 


)P3]} 


(5-2) 


Isothermal  integration  of  — from  the  reference  to  the 

P 

final  density  yields  the  change  in  the  Hemholtz  free  energy 


A ( p , T ) -A  ( p , T ) p p P . . 

= (-^-1)  - In— -2. — (- — — ) + C*{  (9.8642- 


RT 


p*RT  p p 


[ ( 9 . 8 6 4 2 .10^.  K5356  ) . ( : 4 . 2 , 2 5 . 
T T2  p T T2 
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)po+(^*18067  2j_9^43_3  )p^_(2. 06515- 


3.18425  1.00425  ^-3-,  ,po  , x , 

1 Z2 )P0](- — 1 ) + ( 14.2325- 


15.432  3.0147w . 

— Z2 ' (P-P0)  + 


(4.59033--^iLg2-l:.45217)(~2_~2)_(0<68838_lI06142_0_._33475) 

T ?2  ° r T2 


<p3-5^)} 


(5-3) 


Isobaric  temperature  differentiation  of  the  pressure 
yields  the  thermal  pressure  coefficient,  yv . 
dP  P-P  dp 

Y v " ( — °)  + ( -)  - RT( — -)  - p*C*R{-[(10.191+3,0712) 

dT  T dT  T T 

(P-P  Q)- (15. 432+^il29*)  (p2-p2)  + (l0.966+^l°i^)(p3-p3) 

T T 0 

-(3.18425+2l°^)(p~4-p4)]  - T[(9.8642-i£4-191-^5336) 


(-^-2-) -(14.2325 — 15,432 — 3-,0147)(-dg°  ) + 


dT 


dT 


dp3 


(9.18067  966  -— * 9-°  ^ 3 3 ) ( — °)-(2. 06515-  3,18425 


dT 


1^00425)  (^C)]} 

T 2 dT 

(5-4) 

where  the  derivatives  of  the  reference  properties  PQ  and 

P0  are  along  the  path  chosen  for  their  variation,  usually 
along  saturation  line. 
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The  isobaric  expansivity  ap  then  can  be  calculated  through 
the  relation 

aP  = hT^v  (5-5) 

From  the  Maxwell  relation 


/ 3 S v _ / 3 P s 


The  entropy  derivative  is 


(9S) 

jT 


T. 


(5-6) 


(5-7) 


Substitution  of  Equation  (5-4)  into  Equation  (5-7)  and 
isothermal  integration  from  the  reference  density  to  the 
final  density  yields  the  entropy  change. 

S ( p , T ) -S  ( p , T ) i i . d P i i dp 

R - R(p—)(HT“)  “ T(“— )(  — ) + 

o P PQ  dT 

1 - C*{  (9.8642---°-‘382-4,6°68)ln^ (14  . 2?2?-30 ' 864-9  ‘ 0441 


( P - P n ) + ( 4 . 5 9 0 3 3 - ) ( P 2 - P 2 ) - ( 0 . 6 8 8 3 8 - 


1 . 00425  v ,~3  ~3n 


dp. 


L^±^)  (p^-p^-TE  (9.8642-— ■191-1.-  bl56)  (—5.)-  (14.2325- 
T ° T T^  dT 


dp2 


dT 


dT 


^ i 3.18425  1 . 00425  w dpo  ,,,  1 1 , 

(2.06515 ) ( )]( )}  + {l-C*[  (9.8642- 

T T2  dT  p p0 

20,3j2.».«068)_(14 ,2,25-3°-;6A-9^41)-o+(;-180672K922. 


T2 


T 


^4^ ) p 2 - ( 2 . 0 6 5 ! 5 ) ^ ] } ( !£. x ) 
-pZO  ~ “r  Z O''*' 


(5-8) 
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Internal  energy  U is  given  by 


dU  = TdS  - Pdv 


(5-9) 


Isothermal  integration  and  substitution  of  Equations  (5-1) 
and  (5-8)  into  Equation  (5-9)  yield  the  internal  energy 
change . 


U(p,T)-U  (p  T)  Pn  dP  . . dp 


ITT 


o T dT 


P P0  dT 


C^{i(10.191+l^Zi2)ln^_(15>432  + 6JL02^)(~_~  ) + (5< 


483  + 


2-9°433)(p2-^)-(1.06142+£^^)(p3-^)]^[ao.l91+l^Zll) 

-(15.432+11^)5  +(10.966+l^£l«2)p2-(3.1B425+il2£85)53] 

T 0 T 0 T 0 

(— -1)  + Tf  (9.8642-J^jili-ljiiliU-^g_)-(i  /|  2?2?  15’432 
p T T2  dT  T 

2 

3l^)(l^)  + (,.18067.10^_2190«2)(^),(2_Q6515_ 

T dT  T T2  dT 


*+*  q. 

3.18425  1 . 0 0 4 2 5 v,  ,dp0)  -j  (1_1  } j 
T T^  dT  p pQ 


(5-10) 


The  enthalpy  H is  found  from  the  definition 


H = U + PV 


(5-11) 
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Substitution  of  Equations  (5-1)  and  (5-10)  into  Equation 
(5-11)  gives  the  enthalpy, 


H = U0(p0,T)  + 


po(fi->  + T f-rr2) 


P P, 


dP 

( 

'■3T' 


RT[(l-_°)_r(i-i_) 

P P P„ 


dp 


(_£)  + C*RTi-(in.l91^2lLiZlinnP_u.n/i  ?^_30-864  9.0441 


dT 


(P-P  o)-(9.1807-l6--449-l:-8°863)(g2-p2)  + (2.60515- 

T T 0 


4.24567 


1 . 67375  v ,~3  ~3 


tv  ' y > (P^-Pn)  + C (9.8642+illJll)-(i4. 23  25+^-lHiiI)p 
0 y2  Mo 


3.0147 


(9.1807+il££^Z)p2_(2.60515+ll^25)p^](^-l)  + T[(9.8642- 

T T p 


10.191  1.53»«)(^).(14j2325_iL432_3.Q1»7)(<)  + 


dT 


dT 


(9. 18067---0, 966----?-0433  ) Q6515  3,18423  1-00425 


T dT 


dp  . . 

(_£)  ](!-!_)} 

dT  p H 


(5-12) 


In  the  above  equations,  the  property  value  at  the  state  p,T 

can  be  obtained  completely  in  terms  of  the  property  at  the 

state  P0,T  and,  where  necessary,  derivatives  such  as 
dp 

(-j j ) of  that  state. 
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Pure  Liquid  Compression  Data  Bank 
Critical  compilations  of  the  experimental  data  of 
compressed  liquids  are  scarce.  A recent  survey  of  PVT 
properties  of  liquids,  as  referenced  in  Chemical  Abstracts 
before  the  end  of  1983,  was  conducted  by  Tekac  et  al.^ 
and  published  in  the  Journal  of  Fluid  Phase  Equilibria  in 
1985.  This  survey  may  serve  as  a source  of  references  to 
some  original  papers  on  PVT  properties  of  compressed 
liquids.  Our  present  pure  liquid  compression  data  bank 
contains  133  organic  and  inorganic  liquids  (including 
liquefied  gases  and  metallic  liquids),  13  molten  salts,  17 
polymer  melts,  and  7 deuteriated  liquids.  They  are  listed 
separately  in  the  first  column  from  Table  5-1  to  Table 
5-4.  The  substance  order  is  based  on  the  Chemical 
Abstracts  system.  In  the  case  where  the  substance  has  been 
studied  by  several  authors  or  by  the  same  author  with 
several  investigations,  the  reference  numbers  are  ordered 
chronologically.  Only  those  data  sets  with  the  pressure 
ranges  greater  than  300  bar  and  temperatures  more  than  two 
isotherms  have  been  selected.  There  are  234  total  data 
sets.  Several  misprints  have  been  found  in  the  printed 
values  of  the  experimental  data  in  the  original  articles. 
They  have  been  corrected  in  the  present  data  bank. 

Experimental  determination  of  the  volumetric  proper- 
ties of  compressed  liquids  may  be  divided  into  two 
methods:  direct  measurements  and  indirect  measurements. 

The  major  sources  of  the  experimental  data  in  the  present 
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data  bank  come  from  the  direct  measurements  of  volume 
change  or  volume  itself  as  a function  of  pressure  at 
constant  temperature  such  as  from  piezometric^®  or  pycno- 
metric  methods. ^9  Isochoric  data9^  are  generally  not 
included  in  the  data  bank  because  they  are  not  available 
in  isothermal  form.  Only  a small  portion  comes  from  the 
indirect  methods  where  liquid  density  can  be  related  to 
the  other  physical  properties  such  as  ultrasonic  velocity, 
dielectric  constant,  refractive  index,  and  heat  capacity, 
etc.  In  the  case  of  the  ultrasonic  velocity  measurement, 


i Tap 

T = F7  + 


(5-13) 


where  u is  the  velocity  of  sound.  Isothermal  integration 
from  the  reference  pressure  PQ  to  the  final  pressure  P 
yields  the  density  expression.^ 


P = P0  + 


T 

u 


dP  + T 


P 2 
P 


-dP 


(5-14) 


o ro 

In  the  case  of  the  dielectric  constant  measurements,  the 
theoretical  expression  of  the  Clausius- Mossotti  Equation29 
can  be  written  as 


n 1 e-1 

p ■ rn  in 


(5-15) 
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where  CM  is  the  Clausius-Mossotti  function7^  74  and  e is 
the  dielectric  constant.  Similarly,  for  infractive  index 
measurements  the  Lorentz-Lorenz  Equation^  can  be  written 
as 


1 n -1 

p = TT  TT  (5-16) 

n +2 

where  LL  is  the  Lorentz-Lorenz  analogue74  of  the  Clausius- 

Mossotti  function  CM  which  replaces  e in  Equation  (5-15) 

2 

with  n , the  square  of  the  refractive  index. 

The  direct  measurements  tend  to  be  more  accurate  than 
the  indirect  measurements.  The  accuracy  of  the  experimen- 
tal data  varies  from  0.01%  to  1%  in  density.  Generally,  a 
precision  of  0.1%  in  density  corresponds  to  1%  in  pressure 
or  1-10%  in  isothermal  compressibility. 

Results  and  Discussions 

The  pressure  Equation  (5-1)  can  be  used  as  an  equation 
of  state  for  pure  liquids.  If  the  characteristic  parameters 
are  known  for  a specific  substance,  then  the  change  in 
pressure  due  to  an  isothermal  change  in  density  can  be 
calculated.  On  the  other  hand,  Equation  (5-1)  can  also  be 
used  as  a correlation  scheme  to  fit  the  liquid  compression 
data.  Tables  5-1  to  5-4  list  the  characteristic  parameters 
for  all  of  the  substances  in  the  present  data  bank.  The 
results  are  in  excellent  agreement  with  the  experimental 
data  for  all  of  the  substances  over  all  data  ranges.  In 


Characteristic  Parameters  and  Comparisons  with  Liquid  Compression  Data 
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fact,  they  all  appear  to  be  within  experimental  error. 

This  proves  the  generality  and  flexibility  of  the  present 
correlation.  The  characteristic  parameters  for  each  sub- 
stance evaluated  from  the  liquid  compression  data  with 
those  from  the  liquid  DCFI  data  are  very  close  to  each 
other  except  for  the  data  of  Grindley.^"5  This  discrepancy 
is  probably  due  to  the  experimental  uncertainty  of  liquid 
DCFI  at  low  pressure  where  Grindley's  mathematical  treat- 
ment of  the  compression  data  is  questionable.  The  general 
similarity  of  the  characteristic  parameters  for  each  sub- 
stance indicates  the  thermodynamic  consistency  of  the 
present  approach.  Both  the  compressibility  Equation  (4-1) 
and  the  pressure  Equation  (5-1)  with  the  same  parameters 
reproduce  the  crossover  behavior  of  the  volumetric  proper- 
ties of  compressed  liquids  discussed  in  Chapter  3. 

The  characteristic  parameters  are  optimized  only  in 
the  tabulated  data  ranges.  Thus  computations  performed 
outside  the  data  ranges  for  some  substances  may  be  less 
accurate.  This  is  especially  true  for  those  highly  polar 
and  associating  species  where  the  characteristic  parame- 
ters are  evaluated  from  a very  narrow  data  ranges  or  from 
questionable  data  sets  such  as  the  case  of  Schornack  and 

Q C 

Eckert.  J The  characteristic  parameters  are  evaluated 
from  two  isotherms  with  a temperature  span  of  20°K  and  the 
experimental  accuracy  in  density  of  0 . 3 ?o  only.  Though  the 
fitted  results  are  within  the  experimental  errors  (less 
than  3?o  in  pressure),  the  characteristic  parameters  are 
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probably  not  as  reliable  as  those  from  the  more  accurate 
experimental  results.  Even  so,  extrapolations  from  the 
present  correlation  should  be  more  reliable  than  from  other 
correlations.  This  is  particularly  important  in  physical 
properties  correlation  since  experimental  data  covering 
the  whole  state  conditions  are  not  always  available. 

There  are  a number  of  substances  where  several  sets 
of  characteristic  parameters  are  tabulated  together 
because  several  sets  of  measurements  exist.  For  a sub- 
stance with  similar  parameter  sets,  in  most  cases,  any  one 
set  can  be  used  as  the  characteristic  parameter  set.  This 
is  the  case  for  most  of  the  normal  liquids  where  the 
characteristic  parameters  are  found  over  a wide  range  of 
state  conditions.  For  substances  with  different  parameter 
sets,  the  recommended  values  are  characterized  with  an 
asterisk.  Variations  of  the  characteristic  parameters 
are  due  mainly  to  the  precision  of  the  experimental 
results  and  the  insensitivity  of  the  DCFI  and  pressure 
changes  to  parameter  variations.  It  may  be  that  some 
sets  are  not  equivalent  because  of  limited  ranges  of 
measurement  as  for  some  polar  and  associating  substances. 
The  recommended  characteristic  parameters  are  based  on  the 
data  range,  experimental  precision,  and  accuracy  of  the 
measurements.  There  is  a preference  in  this  work  for  the 
data  of  some  authors.  Generally,  they  are  liquefied  noble 
gases  and  others  from  Streett  and  his  coworkers,^’^ 
liquefied  natural  and  petroleum  gases  from  the  National 
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Bureau  of  Standards  (NBS),60  63  hydrocarbons  from  Dymond 
et  a 1. , 1 ^ ’ 1 3 1 ’ 13  3 long  chain  hydrocarbons  from 
Doolittle,138  coal  liquids  (PSU  compounds)  from  Cutler,  et 
and  Lewitz  et  al.,13^  benzene  and  its  derivatives 
from  Gibson  and  Loe f f ler , 47 ’ ^ cyclic  compounds  from  Kuss 
and  Taslimi,  polar  liquids  from  Kumagai  and  his  co- 
workers,78’84’138  alcoholic  liquids  from  Makita  and  his 
coworkers,91’97  liquid  acids  from  Karpela,95  molten  salts 
from  Barton  and  Speedy,144  polymer  melts  from  Simha  and 
his  coworkers,147’149  and  deuteriated  liquids  from  Jonas 
and  his  co w orkers 83 ’ 93 ’ 198-1 94 

Group  Contribution  Method 

In  addition  to  the  molecular  corresponding  states 
correlations  of  the  above  type,  group  contribution  methods 
are  very  useful  for  physical  property  estimation.  The 
concept  assumes  that  each  molecule  is  made  up  of  different 
types  of  atomic  or  structural  groups  and  that  the  molecu- 
lar properties  can  be  evaluated  from  contributions  of 
these  atomic  or  structural  groups.  Unlike  the  correspond- 
ing states  correlation,  the  fundamental  assumption  is  that 
various  groups  in  a molecule  contribute  a definite  value 
to  the  total  properties  independent  of  the  presence  of  the 
other  groups  on  the  molecule.199 

There  are  two  general  approaches  in  group  contribu- 
tion method.  The  first  one  may  be  called  the  property- 
additive  approach.  It  assumes  that  the  groups  actually 
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posses  thermodynamic  properties  or  make  directly  additive 
contributions  to  the  molecular  properties. 15^  The  mole- 
cular properties  then  can  be  calculated  from  the  sum  of 
these  group  properties. 


M.  = I v.  M 
1 ia  a 

a 


(5-17) 


where  v^a  is  the  stoichiometric  coefficient,  represents 
the  number  of  groups  or  type  a in  molecule  i.  This 
approach  has  been  applied  to  the  liquid  molar  volume  at 
the  normal  boiling  point15^  and  to  the  liquid  heat  capa- 
city,  as  well  as  to  the  solution  of  group  methods  for 
activity  coefficient.159’ 160  This  group  contribution  method 
has  also  been  widely  used  for  ideal  gas  properties  such  as 
heat  capacity,  entropy,  and  enthalpy  of  formations155  etc., 
where  intra-  or  intermolecular  forces  play  no  role.  For  the 
DCFI,  Equation  (5-17)  would  yield 


C . 

l 


Z v.  C 
ia  a 
a 


(5-18) 


The  other  general  group  contribution  method  may  be 
called  the  parameter-additive  approach.  It  assumes  that 
the  molecular  property  can  be  written  as  a function  of  the 
state  variables,  x_,  and  a set  of  parameters,  0. 


M.  = f(x;0) 


(5-19) 
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where 


0. 

1 


v 

ia 


a 


0 


a 


(5-20) 


This  approach  has  been  applied  to  the  estimation  of  the 
critical  properties  by  Lydersen  . 

Evaluation  of  the  molecular  parameters  from  Tables 
5-1  to  5-4  shows  that  these  parameters  appear  to  be  group 
additive,  according  to  Equation  (5-20).  Thus,  the  DCFI 
can  be  written  as 


Ci  = C*f(pV*,T*/T) 


(5-21) 


and 


V* 

l 


z 

a 


v . \l* 
la  a 


(5-22) 


C* 

i 


2 

a 


v.  C* 
ia  a 


(5-23) 


Table  5-5  lists  the  group  parameters  that  can  be  added  for 
each  group  to  estimate  the  crossover  parameters  for  some 
of  the  substances.  By  adopting  T*=0.96TC,  comparisons  of 
the  group  contribution  correlation  with  liquid  n-alkane 
compression  data  are  listed  in  Table  5-6.  The  results  are 
satisfactory. 
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Table  5-5 

Group  Contributions  to  Characteristic  Parameters  C*  and  V* 


Group  (Attachment)  -C*(+0.5) 

V*,  cm3/mol(  + 0.05) 

-CH^  (Linear  paraffin)  12.6 

26.83 

-CH^  (nonlinear  paraffin)  9.7 

26.02 

-CH^-  (linear  paraffin)  4.3 

18.35 

— C H 2 — (nonlinear  paraffin)  7.9 

16.54 

-CH^-  (ring)  5.8 

18.39 

-CH<  (nonring)  2.3 

12.20 

-CH<  (aromatic  ring)  6.8 

14.77 

>C<  -3.4 

9.95 

>Si<  (silicon)  11.9 

17.72 

>Sn<  (tin)  16.4 

28.20 

-F  (on  nonring)  6.7 

12.40 

-Cl  (on  nonring)  12.5 

21.29 

-Br  (on  nonring)  9.1 

28.23 

-C00H  (on  nonring)  21.2 

28.64 

-COO-  (ester  on  nonring)  25.0 

21.10 

Group  Contribution  Correlation  and  Comparisons  with  n-Alkane  Liquid  Compression  Data 
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Group  parameters  are  estimated  from  these  data  sets 


CHAPTER  6 

APPLICATION  OF  THE  MODEL  TO  LIQUID  MIXTURES 


Introduction 

Mixtures  of  compressed  liquids  are  the  most  prevalent 
form  of  material  found  in  high  pressure  chemical  pro- 
cesses such  as  hydrogenation,  polymerization,  and  some 
other  energy-related  processes  as  well  as  geological 
systems,  including  oil  and  steam  reservoirs,  for  example, 
the  syntheses  of  ammonia  and  methanol,  petrochemical  or 
coal  hydrogenation,  the  fabrication  of  polyethylene,  etc. 
Despite  their  technological  importance,  the  state  of  know- 
ledge about  their  compositional  and  density  dependence  of 
the  thermodynamic  properties  of  compressed  liquid  mixtures 
is  still  not  known  in  enough  detail  to  be  applicable  to 
process  calculations.  Though  there  exists  a rigorous 
formula  for  the  reduced  bulk  modulus  of  a mixture  in  terms 
of  the  quadratic  mole  fraction  average  of  the  binary 
DCFI,  the  present  model  of  Equation  (4-1)  does  not 
provide  a means  for  obtaining  the  unlike  pair  DCFI.  An 
expedient  alternative  is  to  use  a one-fluid  concept^^^  to 
obtain  the  mixture  characteristic  parameters  as  a function 
of  composition  and  the  pure  component  characteristic  para- 
meters. The  pressure  equation  with  mixing  rules  for  the 
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parameters  has  been  applied  to  all  of  the  compression  data 
in  the  present  compressed  liquid  mixture  data  bank.  The 
results  with  and  without  a binary  parameter  are  tabulated. 
For  mixtures  with  negligible  excess  volume,  the  simple 
(mole  fraction  average)  mixing  rules  work  as  well  for 
mixtures  as  does  the  fitting  of  the  pure  components. ^^2 
This  is  consistent  with  the  findings  of  the  parameter- 
additive  group  contribution  described  in  the  last  chapter. 

Derived  Thermodynamic  Properties 
Thermodynamic  properties  of  liquid  mixtures  are  com- 
pletely defined  by  the  binary  DCFI.  However,  statistical 
mechanics  does  not  provide  a means  to  evaluate  these 
properties  except  for  rigid  body  systems.27  The  experi- 
mental behavior  of  the  like  and  unlike  DCFI  can  be  deter- 
mined through  the  simultaneous  knowledge  of  the  volumetric 
properties  and  activity  coefficient.  From  Campanella,163 

-2 

v^  dlny^ 

1 ' C11  = vhtRT  + x2  dx2 

V2  dln^i 

1 " C12  = V^rT  + xi-d^T 

dlnY2 

vx jRT  + x2  dxj 


T,P 


T,P 


T,P 


(6-1) 


( 6-2a ) 


( 6-2b ) 
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dlnY2 


(6-3) 


VHyRT  + X1  d x 2 T , P 


where  the  volume  variable  of  the  DCFI  is  here  obtained 
from  the  pressure  variable.  These  formal  relationships 
may  cause  some  difficulties  in  the  evaluation  of  the 
binary  DCFI  due  to  the  lack  of  available  good  quality 
experimental  information.  The  alternative  is  to  use  the 
Van  der  Waals  one-fluid  concept  in  the  corresponding 
states  approach.  The  mixture  is  considered  to  be  a hypo- 
thetical pure  liquid  and  the  reduced  DCFI  for  the  mixture 
is  given  as 


( 6-4a ) 


where 


m 


( 6-4b ) 


m 


( 6-4c ) 


T* 

m 


( 6-4d ) 


In  the  case  of  the  binary  parameter  k . ■ = 0 , the  volume 

vJ 

parameter  reduces  to 
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V*  = 2 x.Vf  (6-4e) 

i 

From  Equation  (6-4),  the  bulk  modulus  of  liquid  mixtures 
can  be  written  as 

HTm  = PRTd-Cjp)  (6-5) 

and  the  pressure  expression  for  the  compressed  liquid 
mixtures  is  given  as 


P-P 

irr  = Kn’po)  (6-6) 

where  the  function  g is  the  integral  from  pQ  to  p of  the 
function  f in  Equation  (6-4a).  The  pressure  dependence  of 
the  excess  volume  at  constant  temperature  and  composition  is 
given  as 


3VE  1 v xi 

IP-  " p2(l-Cm)RT  - l p2(l-C.)RT 


(6-7) 


where  the  subscript  i denotes  a pure  component  property. 
This  function  can  reproduce  the  variations  shown  in  Figure 
3-8.  The  partial  molar  volume  can  be  calculated  as 
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- - C™  ) 

vi  " KJW~Jl,PtU.7ii 


7ZV- 

KHV~  T ,_N 


y/wTm(M“)T,V,NjVi 


(6-8) 


where  is  the  number  of  moles  of  component  i. 

For  the  binary  liquid  mixtures,  the  DCFI  can  be 
related  to  the  Helmholtz  free  energy. 


c 1 pV  32A 


'11  “ x,  TFT 


3N: 


t,v,n2 


(6-9) 


r 1 pV  3 2 A 

C22  = - TFT^I 


T,V,N1 


(6-10) 


'12 


p \l  3 A 


’RT  3N13N2 


T,V 


X1  p V 3 2A 
X2  RT  3N2 


1 V 3P 


t,v,n2 


x,  RT  3 N ■ 


t,v,n2 


(6-lla ) 


X_l_  pV  3jVA 
X1  ^ 3N? 


t,v,n1  *1 


1 V 3P 
x7  FT  TFT 


T,V,N1 


(6-llb) 
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Details  of  the  properties  from  the  present  correlation  are 
listed  in  Appendix  B. 

Liquid  Mixture  Compression  Data  Bank 
Experimental  information  on  the  volumetric  properties 
of  compressed  liquid  mixtures  is  important  in  the  funda- 
mental understanding  of  the  statistical  theory,  structure 
effects,  and  intermolecular  interactions  of  solution  at 
high  pressurers.^^  However,  there  are  considerably  fewer 
liquid  mixture  compression  data  than  there  are  for  pure 
components.  The  present  data  bank  contains  34  binary 
systems.  They  are  arranged  in  the  order  of  increasing 
complexity  of  the  intermolecular  interactions  and  listed 
in  the  first  column  of  Table  6-1.  These  include 
simple  liquid  mixtures  from  Streett  and  his  co- 
workers,^’^’^’®®’^  n-alkane  mixtures  from  Dymond  et 
^ ^ ^ mixtures  of  various  Cg  compounds  from  Matsuo 

I I Q 

and  Van  Hook,  mixtures  of  benzene  and  its  derivatives 
from  Takagi  and  his  coworker,115’116  mixtures  of  alcohol 
and  hydrocarbons  from  Ozawa  et  al.,50  and  acqueous  solu- 
tions of  alcohols  from  Makita  and  his  coworkers.^’^  The 
experimental  accuracy  of  the  mixture  data  is  high  except 
for  those  systems  where  the  volumetric  data  are  derived 
from  ultrasonic  velocity  measurements.^^ 

Results  and  Discussions 

The  extension  of  the  pressure  equation  to  mixtures  by 
using  the  mixing  rules  in  Equation  (6-4)  has  been  applied  to 


Binary  Parameters  and  Comparisons  with  Liquid  Mixtures  Compression  Data 
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Table  6-1  (Continued) 
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all  of  the  compressed  liquid  mixtures  in  the  present  data 
bank.  Table  6-1  summarizes  the  results  with  and  without  a 
binary  parameter  k—.  The  predicted  results  from  the  simple 
mixing  rules  are  listed  in  the  fifth  column  of  the  table. 

The  results  are  quite  good.  There  is  a definite  relation- 
ship of  the  prediction  and  the  accuracy  of  correlating  the 
pure  component  compressions.  The  average  absolute  percent 
deviation  (AAE?o)  increases  with  increasing  complexity  of  the 

binary  system. 

For  simple  liquid  mixtures,  the  agreement  between  the 
experiment  and  calculated  pressures  is  as  good  as  the  indi- 
vidual pure  component.  The  Argon-Methane  system  is  the  only 
exception.  The  relatively  large  deviation  is  probably  due 
to  experimental  uncertainty  which  is  not  cited  in  the 
original  papers. ^ There  is  a striking  feature  in  the 
pressure  dependence  of  excess  volume  in  these  binary 
systems.  The  excess  volume  is  negative  at  low  pressures 
and  rises  rapidly  with  increasing  pressure.  At  pressures 
above  a certain  value,  there  is  very  little  variation  of 
the  excess  volume  on  temperature  or  pressure.  The  present 
correlation  does  reproduce  the  behavior  described  above 
though  the  decreasing  value  of  SV^/aP  | J,x  can  be  faster 
or  slower  than  the  data  show. 

For  n-alkane  mixtures,  the  deviations  tend  to  increase 
with  increasing  difference  in  the  chain  length.  Use  of  the 
volume  fraction  instead  of  mole  fraction  average  in  the 
temperature  parameter  improves  the  predicted  results  for  all 
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of  n-alkane  mixtures  tested.  This  is  associated  with  the 
Van  der  Waals  forces  between  the  segments  instead  of  the 
molecules. 16  The  reduced  temperature  of  a binary  mixture 

is  given  as 


T = Z cp.T. 
. ri  1 
1 


where 


(6-12) 


^ = — 2" 


x . V* 

l l 


Z x.V* 

i = l 1 1 


(6-13) 


Prom  Equations  (6-12)  and  (6-13),  the  characteristic 
temperature  of  the  binary  mixture  then  can  be  written  as 


. . TIT2(xl',i+x2',!) 

m “ x V*T*+x  V * T* 
12  12  12 


(6-14) 


The  Principle  of  Congruence^^  works  out  very  well  only  for 
this  type  of  mixture. 

For  mixtures  of  Cg  compounds,  the  present  correlation 
can  indicate  the  presence  of  order  arising  from  the  nature 
of  the  packing  volume  at  compressed  states.  Physically, 
order  corresponds  to  a cohesion  which  decreases  the  volume 

• b # I f n 

of  a liquid.  Evaluation  of  the  deviation  pattern 
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reveals  the  fact  that  the  calculated  pressure  is  systema- 
tically lower  than  the  experimental  pressure  up  to  highest 
pressures  measured.  This  is  because  the  experimental 
volume  is  consistently  lower  than  the  calculated  volume  at 
a given  pressure,  temperature,  and  composition  since  the 
simple  mixing  rules  do  not  include  the  order  that  appears. 
However,  adjustment  of  mixture  volume  by  taking  into 
account  a binary  parameter  in  the  quadratic  mole  fraction 
average  of  the  volume  parameter  improves  the  overall 
deviations.  The  results  are  as  good  as  the  pure  component 
correlation . 

For  nonpolar  mixtures  of  globular  species,  the  pre- 
dicted results  are  comparable  to  those  of  the  pure  com- 
ponent and  are  independent  of  the  size  difference  of  the 
molecules.  This  is  the  case  of  carbon  tetrachloride- 
octamethylcyclotetrasiloxane81  in  which  the  latter  has  a 
characteristic  volume  more  than  three  times  of  the  former. 
For  mixtures  of  benzene  and  its  derivatives,  the  devia- 
tions increase  with  increasing  polarity  of  the  substi- 
tuents. The  large  deviation  is  due  probably  to  the 
experimental  imprecision,  which  is  associated  with  an 
anomaly  observed  in  the  ultrasonic  velocity  measure- 
ments. For  alkane-alkanol  mixtures,  the  predicted 

results  are  acceptable  if  the  nonideality  and  experimental 
inaccuracies  are  allowed  for.  The  experimental  excess 
volume  varies  with  composition  in  a complicated  manner  and 
nonideality  is  significant  only  at  pressures  below  1 kbar. 
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The  deviation  is  proportional  to  the  absolute  value  of  the 
ambient  excess  volume  at  the  given  composition.  For  mix- 
tures of  aniline-nitrobenzene,  the  pure  T*  and  C*  values 
are  very  similar  to  each  other;  the  only  difference 
is  in  V*.  Though  the  predicted  results  have  only  about  1% 
error,  it  is  not  as  good  as  the  accuracy  of  the  experi- 
ental  data.  Use  of  one  binary  parameter  in  the  composi- 
tion dependence  of  volume  parameter  does  improve  the 
result  to  as  good  as  the  pure  component.  For  aqueous 
systems,  the  state  dependence  of  the  volumetric  properties 
is  even  more  complicated  though  predictions  within  10% 
error  are  always  obtained. 

The  results  with  a fitted  binary  parameter  in  the 
composition  dependence  of  the  volume  parameter  are  listed 
in  the  sixth  column  of  Table  6-1.  For  most  of  the  sys- 
tems, the  correlated  results  are  much  improved  in  compari- 
son with  the  predicted  results  and  one  binary  parameter  is 
enough.  For  those  systems  with  complicated  state  depen- 
dence, more  than  one  binary  parameter  would  be  required 
for  agreement  within  experimental  uncertainty.  Though  use 
of  three  binary  parameters  in  the  composition  dependence 
of  the  three  characteristic  parameters  would  completely 
recover  the  same  accuracy  as  the  pure  component  correla- 
tion for  each  composition,  predictive  capability  is 
missing  completely. 


CHAPTER  7 

CONCLUSIONS  AND  RECOMMENDATIONS 


This  work  has  focused  on  the  development  of  a general 
and  accurate  liquid  phase  equation  of  state  which  is  based 
on  an  unique  approach  of  the  statistical  mechanical  direct 
correlation  function  integrals.  Two  treatments  have  been 
adopted  to  model  the  liquid  DCFI:  the  corresponding 

states  principle  and  the  group  contribution  method.  The 
major  engineering  contribution  of  the  present  work  is  the 
construction  of  a three-parameter  corresponding  states 
correlation  which  is  based  on  the  detailed  analysis  of  the 
experimental  DCFI  in  the  compressed  liquid  phase.  The 
molecular  corresponding  states  correlation  either  in  com- 
pressibility form  (Chapter  4)  or  in  pressure  form  (Chapter 
5)  has  been  shown  to  be  the  most  accurate  in  representing 
the  volumetric  properties  of  compressed  liquids  and  liquid 
mixtures . 

In  the  absence  of  a rigorous  site-site  interaction 
model  for  the  direct  correlation  function,  this  work 
developed  an  empirical  parameter-additive  group  contribu- 
tion corresponding  states  correlation  arising  from  the 
appearance  of  the  group-additive  crossover  parameters. 

The  group  contribution  corresponding  states  correlation 
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has  been  shown  to  be  successful  in  the  prediction  of  the 
volumetric  properties  of  compressed  liquid  n-alkanes  in 
which  the  molecular  characteristic  temperature  can  be 
estimated.  However,  application  to  the  other  substances 
is  still  hampered  by  a lack  of  knowledge  about  the 
characteristic  temperature.  In  practical  use,  the  alter- 
native is  to  make  one  compression  measurement  to  determine 
the  characteristic  temperature  while  characteristic  volume 
and  DCFI  can  be  estimated  from  the  group  parameters  listed 
in  Table  5-3. 

An  important  by-product  in  the  process  of  the  group 
contribution  formulation  is  the  establishment  of  an  exten- 
sive data  bank  for  the  volumetric  properties  of  compressed 
liquids  and  liquid  mixtures.  The  experimental  PVT  data, 
either  from  direct  or  indirect  measurements,  have  been 
carefully  evaluated  and  converted  to  common  units,  bar- 
cm^/mol-°k. 

The  most  important  scientific  contribution  of  this 
work  is  the  finding  of  the  crossover  behavior  of  the 
compressed  liquid  isotherms.  For  every  liquid,  there 
exists  a unique  density  where  the  reduced  bulk  modulus  is 
independent  of  temperature.  Thus  the  compressional  pro- 
perties of  a real  liquid  do  not  behave  like  those  of  a 
rigid  body  fluid.  Furthermore,  the  assumption  of  linear 
density  behavior  of  dense  liquid169’170  is  not  fully 
correct.  In  fact,  while  there  exists  a transition  region 
of  density  in  which  the  volumetric  behavior  of  compressed 
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liquid  may  be  linear  in  density,  below  or  above  this 
region  the  corrections  for  the  nonlinear  density  behavior 
have  to  be  introduced. 

Thermodynamic  properties  of  major  interest  have  been 
derived  from  the  present  correlation.  For  pure  liquids, 
the  compressional  thermodynamic  properties  can  be  pre- 
dicted from  the  data  along  the  orthobaric  curve.  For 
liquid  mixtures,  the  excess  thermodynamic  properties  can  be 
calculated  along  with  specific  mixing  rules  for  the 
known  pure  component  properties.  The  correlation  method 
also  provides  a way  to  estimate  the  binary  DCFI  which  can 
be  related  to  the  thermodynamic  derivative  properties  in 
the  fluctuation  solution  theory. 

The  above  conclusions  indicate  our  present  concept  of 
compressed  liquid  properties.  There  are  some  suggestions 
for  future  work. 

1.  To  rework  the  form  of  the  polynomial  to  have  all 
isotherms  pass  through  the  crossover  point. 

Currently,  extrapolation  of  the  present  correlation 
to  the  temperature  below  the  freezing  point  of 
methane  does  vary  away  from  this  point.  The  revised 
formula  will  be  given  as 

m n 

=2  2 aii  TJ_1(p-l)1 

i=l  j=l  J 


C*-C 

TTZ* 


(7-1) 


105 


This  ensures  that  all  isotherms  will  pass  through  the 
crossover  point.  Application  to  low  temperatures 
(TC0.5)  would  be  more  reliable  with  the  modified 
correlation . 

2.  To  investigate  thoroughly  the  pressure  dependence  of 
excess  volume  for  compressed  liquid  mixtures.  An 
empirical  correction  which  could  be  used  to  estimate 
the  excess  volume  of  compressed  liquid  mixtures  from 
that  of  atmospheric  measurement  would  improve  the 
accuracy  of  the  present  correlation,  especially  for 
those  mixtures  with  complicated  variations  of  the 
excess  volume. 

3.  To  apply  the  present  correlation  to  the  calculation 
of  partial  molar  volume  and  chemical  potential  and 
compare  with  experimental  data  for  compressed  liquid 
mixtures.  This  would  be  a strong  test  of  the  model's 
reliability . 

4.  To  optimize  the  present  group  parameters  and  deter- 
mine the  molecular  characteristic  temperature  for 
the  rest  of  substances.  While  the  relative  insensi- 
tivity of  dense  liquid  results  to  the  characteristic 
temperature  means  the  present  group  contribution 
corresponding  states  method  may  find  useful  appli- 
cation in  predictive  work,  refinement  of  the  group 
method  would  ensure  its  value. 


APPENDIX  A 

PERCUS-YEVICK  HARD  SPHERE  DIRECT  CORRELATION  FUNCTION 
INTEGRAL  FROM  VERLET-WE IS  ALGORITHM 

o 

Verlet  and  Weisz  have  presented  a simple  algorithm  to 
calculate  the  effective  hard  sphere  diameter  d(p,T)  from 
Weeks-Chandler-Andersen  (WCA)  perturbation  theory. 

Define  a temperature  dependent  hard  sphere  diameter  as 


dR(T) 


00 


[1  - 


u0(x) 


ITT 


] dx 


( A — 1 ) 


where  x=r/a,  r is  the  separation  distance,  a is  the 
collision  diameter,  and  uQ(x)  is  the  WCA  reference 
potential.  Verlet  and  Weis  show  that  this  can  be 
simulated  by  the  empirical  expression. 


A 0.3837+1.0680 

R O\'4293+0 — 


( A — 2 ) 


where 

0 = e/kT 
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The  density  dependence  of  the  WCA  values  of  hard  sphere 
diameter  to  first  order  is  given  by 

d r dR ( 1 + Y 6 ) ( A — 3 ) 


where 


6 = 71  'O'.  21+404.60 


(A-4) 


Y = 


1-4. 25n  +1.362n2  -0.8751n3 

[to [to [to 


( A-5 ) 


to 


= n - 


I6r 


( A — 6 ) 


it  ,3 

0 = gPd 


( A — 7 ) 


where  n is  the  packing  fraction.  A few  iterations  are 
sufficient  to  find  d.  Then  the  Percus-Ye vick  hard  sphere 
DCFI  can  be  calculated  as 


.pyc  _ n ( 2+n  ) (4-n ) 

'hs  (77? 


( A — 8 ) 


where  the  superscript  PYC  is  the  Percus-Yevick  compressi- 
bility equation  and  the  subscript  hs  is  hard  sphere. 


APPENDIX  B 

THERMODYNAMIC  PROPERTIES  OF  COMPRESSED 
LIQUID  MIXTURES  FROM  DCFI  MODEL 


From 

Equation  (6-4),  the  DCFI  model  can  be  written  as 

C = 

~ ~ ~3 

aQ  + a^p  + a2P  + a^p  ( B — 1 ) 

where  the 

reduced  density  coefficient  a^  are 

ao  = 

9.8642  - 10 . 191x  - 1.5356t2 

II 

f— 1 

CO 

-28.465  + 30.864X  + 6.0294t2 

a2  = 

27.542  - 32.898X  - 8.7130x2 

a3  = 

-8.2606  + 12.737X  + 4.0170x3 

and 

The  pressure  equation  is  given  as 
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p = po  - p*RT{(p-  p0 )-C*[ aQ (p-  P0)+ji(p2- 


5^) 


( B-2 ) 


The  change  of  Hemholtz  free  energy  is  given  as 


A-A 

o 

TT 


p^irr 


i p 

) - C*[aQln— 

P P0  P 


(p-5a) 


^Cp2-t20)^Cp’-el)l 


( B — 3 ) 


For  binary  mixtures,  the  simple  mixing  rules  are  given  as 


r = 

X1VI  + X2V2 

(B-4) 

c*  = 

X1CI  * X2C2 

( B — 5 ) 

T*  = 

X1T1  * x2T2 

( B-6 ) 

Thermodynamic  properties  of  mixtures  then  can  be  obtained 
through  the  following  relationships. 
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where 


3 2 A 


3V‘ 


T,N 


3 P 

"3V 


T,N 


Nh. 


pRT 


[l-C*(a0+a1p+a2p2+a3p3 ) ] 


( B — 7 ) 


3 A / 3 A 


t,v,n2 


T,N 


3 P 
TFT 


pv 


T,V,N2  = "N>T " 


RT 


= -y-1  [l-(C|ao  + C*a^)-2(C*aoP  + C*a1P1  + C*a’p) 


(1+r2-^-^Cia2P2+2C*a2PPl+C*a2p2)  (1+^+-7)  - 

P p p 


i(C*a3p,+3C*a;5p2p1+C*a^p3)(l<_£+-^+_4)  1(1-—) 


5 P2  P3 
o Ko  Ko 


-2  ~3 
P P P 


( B — 8 ) 


a*  = 9.8642  - 10.191t1  - 3.0712tt1 
a[  = -28.463  + 30.864t1  + 12.0588tt;l 
a2  = 27.542  - 32.898x1  - 17.426tt1 
a3  = -8.2606  + 12.737x^  + 8.0340x1^ 


Ill 


and 


"! 

v~ 


T! 


From  Equations  C B — 5 ) and  (B-6),  the  partial  molar  volume 
be  calculated  as 


v - (3V  ) 

1 " ^'3TTJ;T,P,n2 

(9P  ) 

= -Tsrp— 

vy\r;T,N 


Ptl”C*(3Q-f-32P  + 32P  *+■  3 P ) 1 


[l-(C*a0+C*a»)- 


1 

2 


(C*a1p+C*a1p1+C*a ’ p ) (1 


)- 


P 


1 

3 


( r* 

^ 1”2 


p ~2 

a„p2+2C*a2pp1+C*a^p2) (l+_£+_^)_ 


i(Cja3p3  + 3C*a3p2p1  + C*a'p3)(l+^+^4j)]  (B-9) 


can 
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M1  1 3 A 

FT  " RT  3N. 


t,v,n2 


= [l-(CJaQ  + C*a^)]lN£_  - [ 1- ( C*aQ  + C*a ’ ) ■ 

P n 


i(CjaiPoH.C*a1p’+c.aiPo) 


-J(CIa2Po+2C*a25o5o+C*a25o)- 


?(cI*3Po+«*a3pJp>+C*a'p’)]a-!2)- 


2(CJa1p+C*a1p,+C^a;[p)(l-^£)- 


i(C*a2p2  + 2C*a2pp1  + C^a’p2)(l-^)- 

P 


j^(C*a  1p'5  + 3C*a;Jp2p1 +C*a^p3  )(l--y) 


(B-10) 


where 
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3 2A 
3N? 


3^ 

T,V,N, 


T,  V,N. 


= 5v{(2C!ai+C*aS)ln^[(2CIa;+C*aJ)+^(2C*aiPi;+2C*aiPo 


+2Ciai5i+C*alPa)+;C4CIa2Po5^2C*a^2+2c*a2Pi;2  + 


4C*a^op;+C.a2p2)+2(6C*a3p2p'+2Ci[a^^6C.a3Pop;2+ 


6C!ajPoPi+C*a3Po>:l(1‘^)-2(2CIal5l+2C!aiP+c*alP>- 


oN  1 


(l~)-g(4cia2PPi  + 2cfa2P2  + 2C*a2Pl  + 4C*a2^^1  + C*a2^Z) 


P2 

^"^7^  ~Y2-(6C*a^p2p^  + 2C*ajP'5  + 6C*a  jPp2  + 6C*a^p2p1  + C*a^ 


P3 

(1-Tj) 

P 


( B — 1 1 ) 


where 
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2 

a”  = 9.8642  - 10.191t1  - 3.0712t1 
a = -28.465  + 30.864x1+  12.0588x2 
a 2 = 27.542  - 32.898t1  - 17.426T2 
a'j  = -8.2606  + 12.737x2  + 8.0340x2 

Comparison  of  Equations  (B-7),  (B-8),  and  (B-ll)  with  the 
following  equations, 


9P 

~a¥ 


T,N 


^^^.(l-x^C  -2x  x C -x^C  ) 

1 1 11  1 2L12  x2 L22 


( B — 1 2 ) 


a / a a 
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RT/1 

pV[T. 


■Cll> 


( B-14 ) 


The  binary  DCFI  then  can  be  calculated. 


Appendix  C 
Computer  Programs 
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VOLUMETRIC  PROPERTIES  OF  ARGON  USING  THE  STROBRIGE  EQUATION 
REF:  TWU  ET.  AL.  FLUID  PHASE  EQUILIBRIA, 4( 1980) 

IMPLICIT  REAL*8(A-H, 0-Z) 

DIMENSION  A ( 1 6 ) 

DATA  A/1. 31024, -3. 80636, -2. 37238, -0. 798872, 0. 198761, 

1 1.  47014,  -0.  786367,  2.  19465,  5.  75429,  6.  78220,  -9.  94904, 

1 -15.  6162,  86.  6430,  18.  5270,  9.  04755,  8.  68282/ 

DO  2 1=1, 8 
READ  (5,  *>  T 

PM=0. 6890606D0+2078.  76667D0*( (T/83.  80D0>**1.  59817868D0 
1-1.  ODO) 

TR=T / 1 50.  86D0 
TS=1.  2593D0*TR 
WR ITE ( 6, 101 ) T, PM 
101  FORMAT (2X,  2F13.  4/ ) 

TI-1.  ODO/TS 

FI =A ( 1 ) +TI* ( A<  2 ) +TI* ( A ( 3 ) +TI* ( A ( 4 ) +A ( 5 ) *TI*TI ) ) ) 

F2=A(6)+A(7)*TI 

F3=A ( 3 ) 

F4=TI*TI*TI*(A<9)+THMA< 10 ) +A < 1 1 > *TI > ) 

F5=TI*TI*TI*<A( 12)+TI*<A< 13 > +A < 14 ) *TI ) ) 

F6=A( 15)*TI 
DO  1 N=l,  40 
READ  (5,*)  V 
RR  = 1. ODO/ ( 0.  01341D0*V) 

RS=0. 3189D0*RR 

R2=RS*RS 

R3=RS*R2 

R4=RS*R3 

R5=RS*R4 

Dl=-A(  16)*R2 

E1=DEXP ( D1 ) 

F= 1 . ODO+F 1 *RS+F2*R2+F3*R3+F4*R2*E1 +F5*R4*E1 +F6*R 5 
DF=Fl+2.  0D0*F2*RS+3.  0D0*F3*R2+F4*E1*<2.  0D0*RS-2.  ODO 
1*A< 16)*R3)+F5*E1*<4.  0D0*R3-2. ODO*A( 16>*R5)+5.  0D0*F6*R4 
OMC=RS*DF+F 

P=F*RR*13.  41D0*0.  0831434D0*T 
IF( P.  GT.  PM)  GO  TO  1 
C=1 . ODO-OMC 
WR ITF ( 6, 5)  V, P, C 
5 FORMAT (3F 13. 4) 

1 CONTINUE 

2 CONTINUE 
STOP 
END 
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MOLLERUP'3  PROGRAM  TO  GENERATE  THE  DCFI  DATA  OF  METHANE 
PROGRAM  METHERM 
PRINTOUT  ISOTHERMS. 

IMPLICIT  REALMS  (A-H, O-Z) 

COMMON/ONE/AG, AL,  BB,  BE,  DE,  GK,  TTRP, PTRP, DTRP,  TCRT, PCPT, DCRT 
COMMON/THRE/DPDT . D2PDT2,  DLPDLR,  DTSDR,  DTHDR,  DPSIDR,  XB1,  XB2,  XD1,  XD2 
COMMON/SI X/GP i G,  GY,  GX,  GXX 
9 FORMAT  ( IX,  3F13.  4) 

CALL  METHAN 
DO  130  1=1,3 
READ  (5,*)  TT 

PM=0.  117435675+1909.  40*( (TT/90.  68>**1.  85-1.  0) 

IF( TT-TCRT)  103,103,104 

103  DG=DENGAS ( TT  > 

DL=DENLIQ(TT> 

104  DO  120  N=l, 34 
READC5,  *)  W 
DN=1000.  ODO/W 
IF(TT-TCRT)  114,115,115 

114  IF(DN.  GT.  DG.  AND.  DN.  LT.  DL>  GO  TO  120 

115  PP  = DPDRF (TT, DN) 

IF(PP.GT.  PM)  GO  TO  130 
WRITE (6,  9)  W,  PP,GY 

120  CONTINUE 
130  CONTINUE 
STOP 
END 

SUBROUTINE  METHAN 
C NOTE  THAT  PRESSURES  ARE  IN  BARS,  1 ATM  = 1.01325  BAR. 

C ONE  BAR-LITER  = 100  JOULES. 

IMPLICIT  REAL*8  (A-H, O-Z) 

COMMON  Al,  A2,  A3,  A4,  Bl,  B2,  B3,  B4,  Cl,  C2,  C3,  C4,  C5,  C6,  Dl,  D2,  D3,  D4,  D5 
COMMON/ONE/AG,  AL,  BB,  BE,  DE,  GK,  TTRP, PTRP, DTRP,  TCRT, PCRT, DCRT 
COMMON/THRE/DPDT,  D2PDT2, DLPDLR, DTSDR, DTHDR, DPSIDR,  XB1,  XB2,  XD1,  XD2 
COMMON/FIVE/  P,  T,  DEN,  DPDD,  E,  H,  S,  CV,  CP,  W,  WK 
COMMON/EIGHT / DDLDT 
C CONSTANTS  FROM  MARSHGAS  9/3/71. 

WRITE<6,  99) 

99  FORMAT  (31H  •»-**  CH4  ***  FEBRUARY  22  1983  ) 

TTRP=90.  68 
PTRP=0.  117435675 
DTRP=2S.  1472 

C PCRT  = 45. 956467  3AR 

TCRT=1.90.  53 
PCRT  =PSATF ( TCRT ) 

DCRT=10.  15 
AG=.  5 
AL=.  5 
BB=0. 8 
BE=4.  0 
DE=1 . 2 

GK  = 0. 0831434 

CJ  = 100 

WK  = 101325/16.  043 
A 1 =-4.  15452847 
A2=-4.  77390186 
A3=3. 51999044 
A4=4. 25737412 
Sl=l. 74744656 
B2=2.  57679569 
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33=0. 4653S980 
01=5. 10335931 
C2=-12. 42022748 
C3=6. 36124273 
Dl  = l.  38383919 
D2=-3. 54713477 
D3=l. 37408815 
RETURN 
END 

FUNCTION  PSATF < T ) 

C METHANE  VAPOR  PRESSURE  VIA  PRYDZ/GOODWIN  DATA,  NOV. , 1970. 

C NOTE,  PRESSURE  IN  BARS,  1.01325  BAR/ATM. 

IMPLICIT  REAL*8  (A-H, 0-Z) 

COMMON/THRE/DPDT,  D2PDT2,  DLPDLR, DTSDR , DTHDR , DPSIDR,  XB1,  XB2,  XD1,  XD2 
DATA  PTRP  / 0.117435675/,  TTRP/90.  68/,  TCRT/190.  53/,  E/1.  5/ 

DATA  A/4.7732553/,  B/ 1 . 7665879/,  C/-0.  5702812/,  D/1. 3311373/ 

1 XK=1-TTRP/TCRT 
X=( 1-TTRP/T)/XK 
Q=l-X 

IF  (Q ) 2,3,4 

2 PSATF  = 0 
DPDT  = 0 

RETURN 

3 W=0 
W1=0 

GO  TO  5 

4 W = G**E 

W1  = -E*W/Q 

5 DXDT=TTRP/XK/T**2 

Z = X*W 

Z1  = W X*W1 

6 FZ  = A*X  + B*X**2  + C*X**3  + D*Z 

7 FI  = A + 2*B*X  + 3*C*X**2  + D*Z1 

8 PSATF=  PTRP*DEXP ( FZ ) 

DPDT=PSATF*F1*DXDT 

RETURN 

END 

FUNCTION  DENLIG(T) 

C METHANE  SATD. LIQUID  DENSITIES  VIA  GOODWIN. 

C THIS  FUNCTION  REFITTED  VIA  , DENSATLQ,  4/29/71. 

IMPLICIT  REAL*8  (A-H, 0-Z) 

COMMON/EIGHT / DDLDT 

DATA  TCRT/190.  53/,DCRT/10.  15/ 

DATA  A/0.539403/,  B/l.  896635/,  C/-0.  018387/,  E/0.88/ 

1 X = T/TCRT 

Z = 1 - X 

IF ( Z.  LT.  .01)  GO  TO  7 

2 F = -B*X**2/Z 

XP  = DEXP(F) 

3 DFDX  = — E* ( 2 ♦ X/Z)*X/Z 

G = Z**0.  36 

4 G = A*Z  + B*G  + C*XP 
DENLIQ  = DCRT* ( 1+G ) 

5 G1  = C*XP*DFDX  - A - 0.  36*B*G/Z 

6 DDLDT  = DCRT *G1 /TCRT 

RETURN 

7 IF(Z.  LT.  1.  D-09)  Z = l.D-09 
G=Z*-».  36 

G=A*Z  + B*Q 
DENLIQ  = DCRT* ( 1 . +G ) 
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Gl  = -A  - . 36*B*Q/Z 
DDLDT  = DCRT*G1/TCRT 
RETURN 
END 

FUNCTION  TSATF ( DEN ) 

TSAT(DEN)  FORMULA.  METHANE,  CONSTRAINED  TO  TRIPLE  POINT,  IPTS-63. 
YIELDS  ALSO  THE  FIRST  DERIVATIVE  RSP.  TO  RHO>DEN/DTRP. 

EQN.  , (TCRT/T-1 >/AZ  = F ( R ) = U(S)*<1  + A1*LQG(R)  + (R-1)*W(R>), 
WHERE,  U(S)  = ( (S-l )/(ST-l ) >**(8/3),  AZ  = < TCRT/TTRP- 1 > , AND  - 
W(R)  = A2  + A3*R  + . . . + A8*R6. 

IMPLICIT  REALMS  (A-H, O-Z) 

COMMON/ONE/AG,  AL,  BB,  BE,  DE,  GK,  TTRP, PTRP, DTRP,  TCRT, PCRT, DCRT 
COMMON/THRE/DPDT , D2PDT2,  DLPDLR,  DTSDR,  DTHDR,  DPSIDR,  XB1,  XB2,  XD1,  XD2 
DIMENSION  A ( 8 ) 

DATA  A / -0.8142449,  2.4794529,  -6.7598384,  41.05600105, 

1  -133.4449230,  223.4955973,  -181.8267154,  58.8359684/ 

1 AZ  = TCRT/TTRP  - 1 

DSDR  = DTRP/DCRT 
SK  = DSDR  - 1 

2 R=D£N/DTRP 
S=DEN/DCRT 

IF  ( (S-D/SK.  LT.  0.  ) Q=-DABS(  < S-l  > /SIX ) **.  333333 
IF(  (S-l  )/SK.  GE.  0.  ) Q=(  (S-D/SK)**.  333333 
IF ( Q ) 3,9,3 

3 U = G**8 

U1  = 8*DSDR*Q**5/SK/3 

4 V = R -1 
W=0 
W1=0 

DO  6 K=2, 8 

5 W = W ♦ A ( K ) *R** ( K— 2 ) 

Ml  » W1  + ( K-2 ) *A ( K ) *R** ( K— 3 ) 

6 CONTINUE 
G=DLOG(R> 

F = U*(l  + A ( 1 ) *G  + V*W) 

7 FI  = U1  + A ( 1 )* ( U/R  + U1*G ) + U*V*W1  + U*W  + U1*V*W 

8 TSATF=TCRT/( 1+AZ*F> 

DTSDR=- ( AZ/TCRT ) *F 1 *TSATF**2 
RETURN 

9 TSATF  = TCRT 

DTSDR  = 0 
RETURN 
END 

FUNCTION  THETAF (DEN) 

IF  S 1,  U = AG* ( S— 1 ) **3.  IF  S A 1,  U = -AL* ( S-l  ) **3. 

YIELDS  ALSO  THE  FIRST  DERIVATIVE  RSP.  TO  RHODEN /DTRP. 

IMPLICIT  REAL*8  (A-H, O-Z) 

COMMON/ONE/AG. AL, BB,  BE,  DE, GK,  TTRP, PTRP, DTRP,  TCRT, PCRT, DCRT 
COMMON/THRE/DPDT, D2PDT2,  DLPDLR,  DTSDR,  DTHDR,  DPSIDR,  XB1, XB2,  XD1,  XD2 

1 S=DEN/DCRT 

DSDR  = DTRP/DCRT 
IF(S)  9,9,2 

2 Q=S-1 
Q2=Q**2 
Q3=Q**3 

I F ( G ) 3,8,4 

3 U = AG*Q3 

U1  = 3*AG*Q2*DSDR 
GO  TO  5 

4 U = -AL*Q3 
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U1  = -3*AL*G2*DSDR 
5 IF  (U.  LT.  -50.  ) U=-50. 

XP  = DEXP  < U ) 

TS  = TSATF ( DEN ) 

THETAF  = TS*XP 

7 DTHDR  = TS*U1*XP  + DTSDR*XP 
RETURN 

8 THETAF  = TCRT 
DTHDR  = 0 

RETURN 

9 THETAF  = 0 
DTHDR= 1 . E 05 

RETURN 
END 

FUNCTION  PHIF(T) 

XB  = PHI  = X*< l-EXP(-U) >,  U = BB  + BE/X, 

YIELDS  ALSO  DPHI/DR,  DPHI/DX,  AND  D2PHI/DX2. 

IMPLICIT  REAL*8  (A-H, 0-Z) 

COMMON/ONE/AG, AL, BB. BE, DE, GK,  TTKP, PTRP, DTRP,  TCRT,  PCRT,  DCRT 
COMMON/THRE/DPDT, D2PDT2, DLPDLR, DTSDR, DTHDR, DPSIDR,  XB1, XB2, XD1, XD2 

1 X=T /TCRT 
U=BB+BE/X 
Ul=~Bt/X**2 
U2=-2*U1/X 

2 XP=DEXP(-U) 

Z=1-XP 
Z1=U1*XP 

Z2-(U2-U1**2)*XP 

3 PHIF  - X*Z 
XB 1 = X*Z1  > Z 
XB2  = X*Z2  + 2*Z1 

9 RETURN 
END 

FUNCTION  PSIF (T, DEN) 

XD  = PSI  = <1-W*L0G< l+l/W) >/X.  W = (T-TH)/(DE*TCRT>. 

YIELDS  ALSO  XD1  > DPSI/DX,  XD2  > D2PSI/DX2,  AND  DPSI/DR 
IMPLICIT  REAL* 8 (A-H, 0-Z) 

COMMON/ONE/AG, AL,  BB,  BE,  DE,  GK,  TTRP, PTRP, DTRP,  TCRT, PCRT, DCRT 
COMMON/THRE/DPDT,  D2PDT2,  DLPDLR,  DTSDR,  DTHDR,  DPSIDR,  XB1,  XB2,  XD1,  XD2 

1 TH  = THETAF (DEN) 

W = ( T-TH ) /DE/TCRT 

IF < W ) 2,2,3 

2 PSIF=1 
XD1  ~Q 
XD2=0 
DPS IDR=0 

RETURN 

3 X=T /TCRT 
DWDR= -DTHDR /DE/TCRT 
DWDX=1/DE 

4 U=l/X 
DUDX=-U/X 
D2UDX2=-2*DUDX/X 

5 G=DLOG< 1.  ODO+1.  ODO/W) 

WS  = 1+W 
V =•  1-W*G 

6 DVDW  - 1/WS-G 

D2VDW2  = 1/W/WS**2 
PSIF  = U*V 

7 DPSIDR  = U*DVDW*DWDR 
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XD1  = U*DVDW*DWDX  + V*DUDX 

8 XD2  = U*D2VDW2*DWDX**2  + 2*DUDX*DVDW*DWDX  + V*D2UDX2 

9 RETURN 
END 

FUNCTION  DPDRF(T, DEN) 

NOTE,  DPDRF  = PRESSURE,  BAR. 

REDUCED  DERIVATIVE,  DLPDLR  > R*<DP/DR)/P  IS  IN  COMMON. 

<Z-1)*X/R  = F ( T, R ) YIELDS  ( DP/DR ) /DTRP/GK/T  = 1 + ( R2*F1 +2*R*F ) /X 
IMPLICIT  REAL»8  <A-H, O-Z) 

COMMON  Al,  A2,  A3,  A4,  Bl,  B2,  B3,  B4,  Cl,  C2,  C3,  C4,  C5,  C6,  Dl,  D2,  D3,  D4,  D5 
COMMON/ONE/AG, AL,  B3,  BE, DE, GK,  TTKP, PTRP, DTRP,  TCRT, PCRT, DCRT 
COMMON/THRE/DPDT , D2PDT2, DLPDLR, DTSDR, DTHDR, DPSIDR,  XB1, XB2, XD1, XD2 
COMMON/SIX/GP, G, GY, GX, GXX 
COMMQN/FQUF/FRT , DFDX, DFDR1 

1 X=T /TCRT 
S=DEN/DCRT 

DSDR  = DTRP/DCRT 

2 R=DEN/DTRP 
R2=R**2 
R3=R-»*3 

3 XB  = PHIF(T) 

XC  = i/x- 

4 XD  = PS IF (T,  DEN) 

XDD  = DPSIDR 

6 A = Al  + A2*R  + A3*R2  + A4*R3 
AD  = A2  + 2*A3*R  + 3*A4*R2 

7 B = Bl  + B2*R  + B3*R2 
BD  = B2  + 2*B3*R 

8 C = C1*R  + C2*R2  + C3*R3 
CD  = Cl  + 2*C2*R  + 3*C3*R2 

10  U = S-l 
U1  = DSDR 

11  V = Dl  + D2*R  + D3*R2 
VI  - D2  + 2*D3*R 

14  D = U*V 

DD  = U*V1  + U1*V 
F = A + B*XB  + C*XC  + D*XD 

16  DFDR  = AD  + BD*XB  + CD*XC  + D*XDD  + DD*XD 

17  O » 1 + <R2*DFDR  + 2*R*F)/X 
G1=G+1 
GY=1 . ODO-G 
GX=DLOG(S) 

GXX— S -1 
GXS-DLQG ( S ) 

G=1+P*R/X 

DLPDLR  = G/Q 
20  DPDRP  = DEN*GK*T*G 
GP=DPDRF*DLPDLR/DEN 
DFDR 1 -DFDR 
RETURN 
END 

FUNCTION  DENGAS(T) 

C METHANE  VAPOR  DENSITIES  VIA  GOODWIN 

IMPLICIT  REAL*8  (A-H, O-Z) 

DIMENSION  A ( 5 ) 

DATA  A / -7.169797,  -1.705077,  1.461998,  0.416246,  0.524880/ 

DATA  TTRP/90.  68/,  TCRT/190.  53/,  DCRT/10.  15/ 

1 U=  ( TCRT-T  > / ( TCRT-TTRP  > 

W—  CICRT/T-1 >/(TCRT/TTRP-l ) 

3 F = A(1)*W  + A<2)*U**0.  36 
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DO  4 K=3, 5 

4 F = F + A ( K > *U»* < K-2 ) 
9 DENGAS  = DCRT*DEXP<F) 
RETURN 
END 
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NBS  PROGRAM  TO  GENERATE  THE  THERMODYNAMIC  SURFACE  FOR  WATER 
MAIN  PROGRAM 
BLOCK  DATA 

IMPLICIT  REALMS  (A-H, O-Z) 

REAL  P,  Q 

COMMON/ACONST/WM,  GASCON,  TZ,  AA,  ZB,  DZB, YB, UREF, SREF 
COMMON/NCONST/G ( 40 ) , 1 1 ( 40 ) , UJ ( 40 ) , NC 
C0MMGN/ELLC0N/G1,  G2,  GF,  Bl,  B2, BIT, B2T , B1TT,  B2TT 
COMMON/BCONST /P ( 10 ) , Q< 10) 

COMMON/ADDCON/ ATZ ( 4 ) , ADZ ( 4 ) , AAT ( 4 ) , AAD ( 4 ) 

DATA  ATZ/2*64.  Dl,  641.  6D0,  27.  Dl/,  ADZ/3*.  319D0,  1.  55DO/,  AAT/2*2.  D4 
$,  4.  D4,  25.  DO/,  AAD/34.  DO,  4.  Dl,  3.  Dl,  1.  05D3/ 

DATA  WM/ 18.  0152D0/, GASCON/.  461522D0/, TZ/647.  073D0/, AA/1.  DO/, NC/36/ 
DATA  UREF, SREF/-432S.  455039D0, 7.  6180802D0/ 

DATA  Gl,  G2,  GF/11.  DO,  44.  333333333333D0,  3.  5D0/ 

DATA  P/.  7478629,  -.  3540782,  2*0.  , . 007159876,  0.  , -.  003528426,  3*0.  / 

DATA  G/l.  1278334,  0.  , -.  5944001,  -5.  010996,  0.  , . 63684256,  4*0.  / 

DATA  G/-.  53062968 529023D3,  . 22744901424408D4,  . 78779333020687D3 
-. 69030527374994D2,  . 17863832875422D5, -.  39514731563338D5 
*, . 33803884280753D5,  -.  13855050202703D5, -.  25637436613260D6 
*, . 4321257598141 5D6,  -.  34183016969660D6,  . 122231 5641 7448D6 
$, . 1 1 777433655832D7 , -.  217348101 10373D7,  . 10829952168620D7 
25441 998064049D6,  -.  31377774947767D7,  . 5291 1910757704D7 
*, -.  1 3302577 177S77D7,  -.  25109914369001D6,  . 465618261 1 5608D7 
*, -.  72752773275387D7 , . 417742461 48294D6, . 14016358244614D7 
*, -■ 3 1 55523 1 392 1 27D7 , . 47929666384584D7, . 40912664781209D6 
%, -.  1 3626369388386D7 , . 69625220862664D6, -.  10834900096447D7 

$, -. 2272282740 1688D6,  . 33365486000660D6, . 68833257944332D4 
$, . 21757245522644D5, -.  26627944829770D4, -.  7073041 8082074D5 
*,  -.  225D0,  -1.  68D0.  . 055D0,  -93.  ODO/ 

DATA  11/4*0,  4*1,  4*2,  4*3,  4*4,  4*5,  4*6,  4*8,  2*2,  0,  4,  3*2,  4/ 

DATA  JU/2,  3,  5,  7,  2,  3,  5,  7,  2,  3,  5,  7,  2,  3,  5,  7,  2,  3,  5,  7,  2,  3,  5,  7,  2,  3,  5,  7 
2,  3,  5,  7,  l,-3*4,  0,  2,  0,  0/ 

END 

SUBROUTINE  BB (T) 

IMPLICIT  REAL*8  (A-H, O-Z) 

REAL  P, Q 

C0MM0N/ELLC0N/G1 , G2,  GF, 31, B2,  BIT,  B2T,  B ITT , B2TT 
COMMON/ACONST/WM,  GASCON,  TZ, AA,  Z, DZ, Y, UREF, SREF 
COMMON/BCCNST /P ( 10 ) , Q(10> 

DIMENSION  V(10) 

V(  1 )=1. 

DO  2 1=2, 10 
2 V ( I > =V  < 1-1 )*TZ/T 

B 1 =P  < 1 ) +P  < 2 > *DLOG ( 1 . / V ( 2 ) ) 

B2=G( 1 ) 

B1T=P(2)*V(2)/TZ 

B2T=0. 

B1TT=0. 

B2T1  =0. 

DO  4 1=3,  10 
B1=B1+P< I )*V< 1-1 ) 

B2=B2+Q( I )*V( 1-1 ) 

S 1T=B IT- < 1-2  >*P ( I ) *V ( I — 1 )/T 
B2T=B2T-< I-2)*Q< I )*V( 1-1 >/T 
B 1 TT=B  1 TT -*-P  ( I ) * ( I -2 ) **2*V  ( I - 1 ) /T/T 
4 B2TT=B2TT+Q( I ) * ( I -2 ) **2*V ( I - 1 ) /T/T 
B 1 TT  =-R  1 TT-B 1 T/T 
B2TT=B2TT-B2T/T 
RETURN 
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END 

SUBROUTINE  BASE(D,T) 

IMPLICIT  REAL*8  (A-H, O-Z ) 

COMMQN/ELLCON/G1 , G2,  GF,  Bl,  B2,  BIT,  B2T , B ITT,  B2TT 
COMMQN/BASEF/AB, GB, SB, UB, HE, CVB, DPDTB 
COMMQN/ACONST/WM,  GASCON,  TZ,  A,  Z, DZ, Y, UREF, SREF 
Y=. 25*B1*D 
X=1 . -Y 

ZO=(l.  +G1*Y+G2*Y*Y> /X**3 
Z = Z0>-4.  *Y*(B2/B1-GF> 

DZ0=(Gl+2.  *G2*Y)  /X**3+3.  *(  1.  +G1*Y+G2*Y*Y > /X**4 
DZ=DZ0+4.  * < B2/B 1 -GF ) 

AB=-DLQG  < X ) - ( G2- 1 . >/X+28.  16666667D0/ X/X+4. *Y* ( B2/B 1-GF ) 
*+15. 166666667D0+DL0G ( D*T*4.  55483D0) 

GB=AB  *-Z 

BB2TT=T*T*B2TT 

UB=-T*B IT* ( Z- 1 . -D*B2 ) /B 1 -D*T*B2T 
HB=Z+UB 

CVB=2.  *UB+ ( ZO-i . )*( (T*B1T/B1 )**2-T*T*BlTT/Bl ) 

S-D*  < BB2TT-GF*B 1 TT*T*T ) - ( T*B 1 T/B 1 > **2*Y*DZ0 
DPDTB=Z/T+D*(DZ*BlT/4.  +B2T-B2/S1*B1T) 

SB=UB  -AB 

RETURN 

END 

SUBROUTINE  GG(T, D> 

IMPLICIT  REAL*8  (A-H, O-Z ) 

COMMQN/RESF/AR, GR,  SR,  UR,  HR, CVR,  DPDTR 
COMMON/QQQQ/Q, G5 

DIMENSION  QR(ll), QT(IO), GZR(9), QZT(9> 

EQUIVALENCE  ( QR (3) , QZR ( 1 ) ) , ( GT ( 2) , QZT ( 1 ) > 
COMMON/NCONST/G ( 40 ) , 1 1 ( 40 ) , UU ( 40 ) , N 
COMMON/ACONST/WM, GASCON, TZ, AA, Z, DZ, Y, UREF, SREF 
COMMON/ADDCON/ATZ ( 4 ) , ADZ ( 4 ) , AAT ( 4 ) , AAD ( 4 ) 

RT=GASCON*T 
GR< 1 >=0. 

Q5=0. 

G=0.  DO 
AR=0. DO 
DADT=0. 

CVR=0. 

DPDTR=0. 

E=DEXP(-AA*D) 

Q10=D*D*-E 
Q20=l.  DO-E 
GR  ( 2 ) =Q10 
V=TZ/T 
GT( 1 >=T/TZ 
DO  4 1=2,  10 
GR  ( I + 1 ) =QR ( I ) *Q20 
4 QT( I )=QT( 1-1 >*V 
DO  10  1 = 1,  N 
K.=  I I ( I ) + l 
L=JJ( I ) 

ZZ=K 

GP=G( I )*AA*QZR(K— 1 >*QZT(L> 

Q=Q+QP 

Q5=G5+AA*(2.  /D-AA*( 1.  -E* ( K-l ) /Q20 ) > *GP 
AR=AR+G( I >*QZR(K)*QZT(L> /Q10/ZZ/RT 
DFDT=G20**K*  ( 1-L ) -*QZT  ( L+l ) /TZ/K. 

D2F=L*DFDT 
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DPT=DFDT*Q10*AA*K/Q20 
DADT=DADT+G( 1 ) *DFDT 
DPDTR=DPDTR+G( I )*DPT 

10  CVR=CVR+G( I )*D2F/GASC0N 
QP=0. 

Q2A=0. 

DO  20  J=37» 40 

IF(G(J).  EQ. 0.  DO)  GO  TO  20 

K=I I ( J) 

KM= J J ( J ) 

DDZ=ADZ ( J-36) 

DEL=D/DDZ-1. 

IF (DABS (DEL) . LT.  1. D-10)  DEL=1 . D-10 
DD=DEL*DEL 

EX  1 = -AAD ( J-36 ) *DEL**K 
DEX=DE XP ( EX  1 > *DEL**KM 
ATT =AAT ( J-36  > 

TX=ATZ ( J— 36) 

TAU=T/TX— 1. 

EX2=-ATT*TAU*TAU 
TEX=DE  XP ( EX2 ) 

Q10=DEX*TEX 

GM=KM/DEL-K*AAD ( J-36 ) *DEL** ( K-l ) 

FC  T=GM*D**2*Q 1 0 /DD  Z 

Q5T=FCT* ( 2.  /D+GM/DDZ ) - ( D/DDZ ) **2*Q1 0* ( KM/DEL/DEL+ 
( K- 1 > *AAD ( J-36 ) *DEL** ( K-2 ) > 

Q5=G5+G5T*G( J) 

GP=QP  i-G  ( J ) -»FCT 

DADT=DADT-2.  *G(  J) *ATT*TAU*Q 1 0/TX 
DPDTR=-DPDTR-2.  *G  ( J > *ATT*TAU*FCT/TX 
G2A=Q2A+T*G( J)*(4.  *ATT*EX2+2.  *ATT)*Q10/TX/TX 
AR=AR  >-Q10*G(  J)/RT 
20  CONTINUE 

SR=-DADT /GASCON 
UR=AR  »-SR 

CVR=CVR+G2A/GASC0N 

G=Q+GP 

RETURN 

END 

SUBROUTINE  DFIND ( DOUT,  P , D,  T,  DPD ) 

IMPLICIT  REAL*8  (A-H, 0-Z) 

COMMON/QQGQ/QO,  Q5 

COMMON/ACONST/WM,  GASCON,  TZ,  AA,  Z, DZ, Y, UREF, SREF 
DD=D 

RT=GASCON*T 
IF  ( DD.  LE.  0.  ) DD=1.  D-8 
IF ( DD.  GT.  1.  9)  DD=1 . 9 
L=0 

9 L=L+1 

11  IF(DD.  LE.  0)  DD=1.  D-8 
IF ( DD.  GT.  1.  9)  DD=1.  9 
CALL  BASE(DD,  T> 

CALL  QQ(T, DD) 

PP=RT*DD*Z+QO 
DPD=RT* ( Z+Y*DZ ) +Q5 
IF(DPD.  GT.  0.  DO)  GO  TO  13 
IF(D.  GE.  . 2967D0)  DD=DD*1.02D0 
IF(D.  LT.  . 2967D0)  DD=DD*. 98D0 
IF ( L.  LE.  10)  GO  TO  9 

13  DPDX=DPD*1. 1 
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IF ( DPDX.  LT.  . 1 > DPDX.  = . 1 
DP=DABS( 1. -PP/P ) 

IF  (DP.  LT.  1.  D-8)  GO  TO  20 
IF(D.  GT.  .3  . AND.  DP.LT.  l.D-7)  GO  TO  20 
IF  ( D.  GT.  . 7 . AND.  DP.  LT.  1.  D-6)  GO  TO  20 
X= ( P-PP ) /DPDX 

IF ( DABS ( X ) . GT.  . 1 ) X=X*.  1 /DABS ( X > 

DD=DD  !-X 

IF(  DD.  LE.  0.  > DD=1 . D-8 

19  IF( L.  LE.  30)  GO  TO  9 

20  CONTINUE 
DOUT=DD 
RETURN 
END 

SUBROUTINE  THERM ( D,  T ) 

IMPLICIT  REAL*8  (A-H. 0-Z  ) 

COMMON/ACONST/WM,  GASCON,  TZ, AA,  ZB,  DZB, Y, UREF, SREF 
COMMON/QGQQ/QP, QDP 

COMMON/BASEF/AB,  GB,  SB,  UB,  HB,  CVB,  DPDTB 
COMMON/RESF/AR, GR, SR, UR, HR,  CVR, DPDTR 
COMMON/ IDF/A  I,  GI,  SI,  UI,  HI,  CVI,  CPI 

COMMON/FCTS/AD, GD,  SD, UD, HD, CVD, CPD,  DPDT , DVDT,  CJTT,  CJTH 
CALL  IDEAL (T) 

RT=GASCON*T 

Z=ZB+GP/RT/D 

DPDD=RT* ( ZB+Y*DZB  > +QDP 

AD=AB  +AR+A I -UREF /T +SREF 

GD=AD+Z 

UD=UB+UR+UI-UREF/T 

DPDT=RT*D*DPDTB+DPDTR 

CVD-CVB+CVR+CVI 

CPD=CVD+T*DPDT**2/  ( D*D*DPDD*GASCON ) 

HD=UD  ►Z 

SD=SB  >-SR+SI-SREF 
DVDT=DPDT /DPDD/D/D 
CJTT-1.  /D-T*DVDT 
C JTH=-C JTT /CPD/GASCON 
RETURN 
END 

FUNCTION  PS ( T ) 

IMPLICIT  REAL*8  (A-H, O-Z) 

DIMENSION  A ( 8 ) 

DATA  A/-7.  8889166D0,  2.  5514255D0,  -6.  716169D0 
$, 33.  239495D0,  -105.  38479D0,  174.  35319D0, -148.  39348D0 
S, 48. 631602D0/ 

IF(T.  GT.  314.  DO)  GO  TO  2 

PL=6. 35731 18D0-8858. 843D0/T+607.  56335D0*T** ( - 6) 

PS=. 1*DEXP(PL) 

RETURN 

2 V=*T/647.  25D0 
W=DABS( 1. DO-V) 

B=0. DO 
DO  4 1=1, 3 
Z = I 

4 S=B+A( I )*W**( (Z+l.  )/2.  ) 

Q=B/V 

PS=22.  093D0*DEXP(Q) 

RETURN 

END 

FUNCTION  TSAT(P) 
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REAL  itS  P,  PS,  TO,  TSAT,  TDPSDT 
TSAT =0. 

IF(P.  GT.  22.  05)  RETURN 
ft=0 

PL=2.  302585*DL0G(P> 

T <3=372.  83+PL*  ( 27.  7589+PL*(2.  3819+PL*(.  24834+PL*.  0193855)  ) ) 

1 IF(T<3.  LT.  273.  15)  TG=273.  15 
IF  ( TG.  GT.  647.  ) TG=647. 

IF  ( ft.  LT.  8)  GO  TO  2 
WRITE(6, 3)  ft,  P,  PP,  TG 
3 FORMAT ( ) 

GO  TO  8 

2 ft=ft+l 
PP=PS ( TG ) 

DP=TDPSDT  (TG*> 

IF(DABS( 1.  -PP/P).  LT.  . 00001 ) GO  TO  8 
TG=TG* ( 1 . + ( P-PP  > /DP ) 

GO  TO  1 

3 TSAT=1 G 
RETURN 
END 

FUNCTION  TDPSDT (T) 

IMPLICIT  REAL*8  (A-H, 0-Z) 

DIMENSION  A ( 8 ) 

DATA  A/-7. 8S89166D0,  2.  5514255D0, -6.  716169D0 
33. 239495D0, -105.  38479D0,  174.  35319D0, -148. 39348D0 
*, 48. 631602D0/ 

V=T/647.  25 
W=l.  -V 
B=0. 

C=0. 

DO  4 1 = 1,  8 
Z=I 

Y=A(I)*W**< (Z+l.  )/2.  ) 

C=C+Y/W*<.  5-.  5*Z— 1.  /V) 

4 B=BfY 
Q=B/V 

TDPSDT=22. 093*DEXP(Q)*C 

RETURN 

END 

SUBROUTINE  IDEAL(T) 

IMPLICIT  REAL*8  (A-H,  0-Z) 

COMMON/ 1 DF/AI,  GI,  SI,  UI,  HI,  CVI,  CPI 
DIMENSION  C ( 18 ) 

DATA  C/.  1 9730271 018D2,  . 20966268 1977D2, -.  483429455355D0 
i,  . 605743189245D1, 22.  56023885D0, -9.  87532442D0, -.  43135538513D1 
*, . 4581 55781D0,  -.  47754901883D-1 , . 41238460633D-2,  27929052852D-3 

$,  . 14481695261D-4, -.  56473658748D-6, . 16200446D-7, -.  3303822796D-9 
*, . 451916067368D-1 1 , -.  370734122708D-13,  . 13754606823SD- 1 5/ 

TT=T/1.  D2 
TL=DLOG ( TT ) 

GI=-(C(I)/TT+C(2) ) *TL 

HI  = (C(2)+C(  1 )*(  1.  DO-TD/TT) 

CP  I =C ( 2 ) — C ( 1 ) /TT 

DO  8 1=3,  18 

GI=GI  -C< I )*TT**( 1-6) 

HI=HI+C(  I )•»(  I-6)*TT**(  1-6) 

8 CPI=CPI+C ( I >*< I-6)*( I-5)*TT**( 1-6) 

AI=GI  -1 
U I =H I -1 
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CVI=-CPI-1 
SI=UI  -AI 
RETURN 
END 

SUBROUTINE  SECVIR ( T, VIR ) 

IMPLICIT  REALMS  (A-H, O-Z) 

COMMON/NCONST/G ( 40 ) , 1 1 ( 40 ) , J J ( 40  > , NC 
C0MMQN/ELLCQN/G1 , G2,  GF,  BB1,  BB2,  BIT,  B2T,  B ITT,  B2TT 
CQMMON/QQGG/QO,  Q5 

COMMON/ACONST/WM,  GASCON,  TC,  AA,  Z, DZ, Y, UREF, SREF 
CALL  BB(T> 

V=TC/T 
VIR-BB2 
DO  3 J=l,  NC 

IF ( 1 1 ( J) . NE.  0 ) GO  TO  3 
L=UU(U) 

VIR=VIR+G( J)*V**(L-1 ) /T/GASCON 
3 CONTINUE 
RETURN 
END 

SUBROUTINE  CORR  ( T,  P,  DL,  DV,  DELG  > 

IMPLICIT  REAL*8  (A-H, O-Z) 

COMMON/QQQG/QOO, Q1 1 

COMMON/ACONST/WM, GASCON, TZ, AA,  ZB, DZB, YB, UREF, SREF 
COMMON/FCTS/AD,  GD,  SD,  UD,  HD,  CVD,  CPD,  DPDT, DVDT, CJTT, C J 
DLIG=DL 

IF( DL.  LE.  0.  ) DLIQ=1.  11-.  0004*T 
CALL  3B(T> 

RT=GASCON*T 

CALL  DFIND(DL, P, DLIQ, T, DQ) 

CALL  THERM (DL,  T) 

GL=GD 

DVAP=-DV 

IF(DV.  LE  0.  ) DVAP=P/GASCON/T 
CALL  DFIND(DV, P, DVAP, T, DQ> 

IF  ( DV.  LT.  5.  D-7)  DV=5.  D-7 
CALL  THERM (DV,  T) 

GV=GD 

DELG=GL-GV 

RETURN 

END 

SUBROUTINE  PCORR ( T, P , DL,  DV ) 

IMPLICIT  REAL*8  (A-H, O-Z) 

COMMON/ACONST/WM, GASCON,  TZ, AA,  ZB,  DZB, YB, UREF, SREF 
P=PS ( T > 

2 CALL  CORR ( T,  P,  DL,  DV,  DELG) 

DP=DELG*GASCON*T / ( 1 . /DV-1.  /DL) 

P=P+DP 

IF( DABS (DELG ) . LT.  1.  D-4)  RETURN 

GO  TO  2 

END 

SUBROUTINE  TCORR ( T,  P , DL,  DV ) 

IMPLICIT  REAL*8  (A-H, O-Z) 

COMMON/ACONST/WM,  GASCON,  TZ,  AA,  ZB, DZB, YB, UREF, SREF 
T=TSAT ( P ) 

2 CALL  CORR (T,  P,  DL,  DV,  DELG) 

DP=DELG*GASCON*T/( 1.  /DV-1.  /DL) 

T=T*(1.  -DP/TDPSDT (T) ) 

IF (DABS (DELG) . LT.  1.  D-4)  RETURN 
GO  TO  2 
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END 

SUBROUTINE  UNIT 
IMPLICIT  REAL*8  (A-H, O-Z) 

DOUBLE  PRECISION  NT, ND, NP, NH, NNT, NND, NNP, NNH 
COMMON/UNITS/ IT,  ID,  IP,  IH, NT, ND, NP, NH, FT, FD, FP, FH 
DIMENSION  FFD ( 4 ) , FFP(5>,  FFH ( 6 ) , NNT ( 4 ) , NND ( 4 ) , NNP(5), NNH ( 6 ) 

DATA  FFD/ 1 . D-3,  1.  DO,  . 0180152D0,  . G16018D0/  , 

DATA  FFP/1.  DO,  10.  DO,  9.  869232667D0,  145.  038D0,  10.  1971D0/ 

DATA  FFH/2*1.  DO,  18.  0152D0,  . 23884590D0, 4.  30285666D0,  . 42D0/ 

DATA  NNT/1HK,  1HC,  1HR,  1HF/ 

DATA  NND/6HKG/M3  , 6HG/CM3  , 6HMQL/L  , 6HLB/FT3/ 

DATA  NNP/6H  MPA  , 6H  BAR  , 6H  ATM  , 6H  PSI  , SHKG/CM2/ 

DATA  NNH/ 6HKJ/KG  , 6H  J/G  , 6HJ/M0L  , 6HCAL/G  , 7HCAL/M0L, oHBTU/LB/ 
DATA  Al,  A2,  A3,  A4/8HTEMPERAT , 7HDENSITY 
$, 8HPRESSURE,  8HENERGY  / 

WRITE<6, 11 > Al 
30  WR ITE ( 6, 12) 

READ (5,  *, END=99)  IT 


IF ( IT.  EQ.  0)  STOP 
IF ( IT.  GT.  4)  GO  TO  30 
NT=NNT ( IT) 

WRITE (6,  11 ) A2 
31  WRITE<6,  13) 


READ (5,  *,  END=99)  ID 

IF ( ID.  GT.  4 . OR.  ID.  LT.  1 ) GO  TO  31 

ND=NND( ID) 

FD=FFD< ID) 

WRITE (6,  11)  A3 

32  WRITE (6, 14) 

READ <5, *, END=99>  IP 

IF<  IP.  GT.  5 . OR.  IP.  LT.  1 ) GO  TO  32 

NP=NNP( IP) 

FP=FFP( IP) 

WRITE (6, 11 ) A4 

33  WRITEC6,  15) 


READ <5, *,  END=99>  IH 

IF C IH.  GT.  6 . OR.  IH.  LT.  1 ) 

NH=NNM( IH) 

FH=FFH( IH) 

RETURN 
99  STOP 

11  FORMAT ( 

12  FORMAT ( 

13  FORMAT ( 

14  FORMAT ( 

15  FORMAT ( 

*ES/MOL, 

END 

FUNCTION 


ENTER  UNITS 
CHOOSE  FROM 
CHOOSE  FROM 
CHOOSE  FROM 
CHOOSE  FROM 
6=3TU/LB  ' ) 

TTT(T) 


GO  TO  33 


CHOSEN  FOR  ' , A8  ) 

1=DEG  K,  2=DEG  C,  3=DEG  R,  4=DEG  F') 
1=KG/M3,  2=G/CM3,  3=M0L/L,  4=LB/FT3') 
1=MPA,  2=BAR,  3=ATM,  4=PSIA,  5=KG/CM2 ‘ 
1=KU/KG,  2=U/G,  3=J/M0L,  4=CAL0RI 


DOUBLE  PRECISION  T,  TTT,  FT,  FD,  FP,  FH,  NT,  ND,  NP,  NH 
COMMQN/UNITS/ IT , ID,  IP,  IH, NT, ND, NP, NH, FT, FD, FP,  FH 
GO  TO  <1,  2,  3,  4),  IT 


1 TTT— T 


FT=1. 


RETURN 


2 TTT =T +273.  15D0 
FT=1. 


RETURN 

3 TTT =T/ 1 . 8D0 

FT=.  55555555 5 5556D0 
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RETURN 

4 TTT=(T+459.  67D0>/1.  8D0 
FT=.  5555555555556D0 
RETURN 

END 

FUNCTION  TTI ( T > 

DOUBLE  PRECISION  T,  TTI,  FT,  FD,  FP,  FH,  NT,  ND,  NP,  NH 
COMMON/UNITS/IT,  ID,  IP,  IH, NT, ND, NP, NH, FT, FD, FP, FH 
GO  TO  (5,  6,  7,  8),  IT 

5 TTI~T 
RETURN 

6 TTI =T-273.  1 5DO 
RETURN 

7 TTI=T*1.  SDO 
RETURN 

8 TTI=T*1.  8D0-459.  67D0 
RETURN 

END 

C MAIN<H20> 

IMPLICIT  REAL*8  (A-H, O-Z) 

DOUBLE  PRECISION  NT,  ND, NP, NH 

COMMON/UNITS/IT, ID, IP, IH, NT, ND, NP, NH, FT, FD, FP, FH 
CQMMGN/QGQQ/QO, Q5 

COMMON/FCTS/AD,  GD,  SD,  UD,  HD,  CVD, CPD, DPDT 
COMMON/ACONST/WM, GASCON, TC, AA, Z, DZ, Y, UREF, SREF 
COMMON/NCONST/G ( 40 ) , 1 1 ( 40 ) , JU ( 40 ) , NC 
DATA  NS1,  NS2/2H  M,  2HFT / 

CALL  UNIT 
NS=NS1 

IF ( ID.  EQ.  4)  NS=NS2 
15  READ (5,  *, END=9)  IOPT,  XISO 
IF ( IOPT.  EG.  0)  GO  TO  9 
GO  TO  ( 101, 201, 301 ),  IOPT 
101  READ (5,  ■*,  END=9)  JOPT,  Yl,  Y2,  YI 
IF(JOPT-l)  15,101,103 
103  TT=X ISO 
T=TTT(TT) 

IF( JOPT.  EQ.  2)  DGSS=Y1 /FP/T/.  4 
IZ=0 

CALL  BB  ( T > 

PSS= 20000. 

DW=0. 

IF(T.  IT.  TC>  CALL  PCORR  ( T,  PSS,  DLL,  DW  ) 

DGSS=DW 

IF  ( DGSS.  EQ.  0.  > DGSS=1.  11-.  0004*T 
PSAT=PSS*FP 

IF(JOPT.  EG.  2 .AND.  Yl.GT.  PSAT)  IZ=3 

IF(Y1. GT.  PSAT)  DGSS=DLL 

IF( JOPT.  GE. 3)  IZ=3 

PIN=Y1-YI 

DIN=P IN 

PINC=YI/FP 

DINC=YI*FD 

22  IF< JOPT.  EQ.  2)  PIN=PIN+YI 
IF(  JOPT.  EG.  3)  DIN=DIN-t-YI 
IF  ( JOPT.  EQ.  2 .AND.  PIN.  GT.  Y2)  GO  TO  101 
IF(JOPT.  EQ.  3 .AND.  DIN.  GT.  Y2)  GO  TO  101 
IF( JOPT.  EQ. 2)  PRES=P IN/FP 
IF(  JOPT.  EQ.  3)  D=DIN/FD 
24  CONTINUE 
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26 

27 


IF (JOPT. EQ. 3 

. OR. 

(JOPT 

IF  (JOPT.  EQ.  1 
TSAVE=  TT-Y I 

. AND. 

T.  LT 

IF  (JOPT.  EQ.  1 
PSAVE=P IN-YI 

. AND. 

IOPT 

IF( JOPT.  EQ.  2 

. AND. 

PRES 

IF(  JOPT.  EQ.  1 

. AND. 

T.  GT 

2 . AND.  PIN. 
GO  TO  26 


LT.  PSAT) > GO  TO  26 


23 


EQ.  2 .AND.  IZ.LE.  2)  TT=TT  I ( TS ) 

GT. PSAT/FP  .AND.  IZ.GE.  2)  GO  TO  26 
TS  .AND.  IZ.GE.  2)  GO  TO  26 


PRES=PSAT/FP 

T=TS 

JOPT.  EQ.  1 ) 
JOPT.  EQ.  2) 
JOPT.  EQ.  2) 
JOPT.  EQ.  1 ) 


IZ=IZ>-1 
IF  (JOPT.  EQ.  2) 

IF(  JOPT.  EQ.  1 > 

IF ( I Z.  EQ.  1 . AND. 

IF(  IZ.  EQ.  1 . AND. 

IFC IZ.  EQ.  2 . AND. 

IF ( IZ.  EQ.  2 . AND. 

CALL  BB(T> 

I F ( IQPT.  NE.  3 
CALL  GQ  ( T,  D ) 

CALL  BASE(D, T> 

RT=GASCON*T 
PDUM=D*RT*Z+QO 
IF  (JOPT.  EQ.  3 .OR.  JOPT.  EQ.  3) 
IF(ICPT.  EQ.  3 .OR.  JOPT.  EQ.  3) 
DGSS=D+PINC/DQ 
CALL  THERM (D,  T) 
U=UD*T*GASCON*FH 
C=DSQRT ( DABS ( CPD*DQ*1 . D3/CVD ) ) 
IF ( ID.  EQ.  4 > C=C*3.  280833 
H=HD*T*GASCON*FH 


DGSS=DLL 

DGSS=DW 

DGSS=DLL 

DGSS=DW 


AND.  JOPT. NE. 3)  CALL  DFIND(D, PRES,  DGSS,  T, DQ) 


PRES=PDUM 

DQ=RT*(Z+Y*DZ)+Q5 


S=SD*GASCON*FH*FT 
DPDTX=DPDT*FP*FT 
DPDD=DQ*FP*FD 
C0MP=1.  D3/D/DQ/FP 
DDDTL-1 . D3*DPDT/D/DQ 
CP=CPD*GASCON*FH*FT 
CV=CVD*GASCQN*FH*FT 
VL=FD/D 

DQUT=1.  0D03*VL 
POUT  =PRES*FP 


ROUT=DOUT / 1 7.  857142D0 
X0UT=R0UT-1.  DO 
GP=DPDD/0.  0831434D0/TT 
GP 1 =1 . DO-GP 

WRITE(6, 21)  TT, POUT, DOUT, GP 1 
21  F0RMAT(F9.  3,  3F12.  5) 

IF(IZ.  EQ.  1)  WRITE(6»  12) 

12  FORMAT ( ' 

* ) 

IF(  IZ.  EQ.  1 ) GO  TO  23 

IF(IZ.  EQ.  2 .AND.  JOPT.  EQ.  2)  PIN=PSAVE 
IFdZ.EQ.  2 .AND.  JOPT.  EQ.  1 ) TT=TSAVE 
IFdZ.EQ.  2)  I Z=3 
GO  TO  (22, 204, 304), IOPT 

201  J0P1=1 
PRES=  X ISO/FP 

202  READ (5,  *,  END=9)  T1,T2,  YI 
IF(T1.  EQ.  0.  ) GO  TO  15 
TT=T1  -YI 

T=TTT ( TT ) 

CALL  TCORR(TS,  PRES,  DLL,  DW) 

D=DLL 
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IF< <T+YI*FT>.  LT.  TS ) IZ=0 
204  TT=TT+YI 
T =TTT ( TT ) 

IF  ( TT.  GT.  T2  > GO  TO  202 
CALL  BB  < T > 

DGSS-D 
GO  TO  24 

301  J0PT=1 
D=X ISO 

302  READ <5<  *,  END=9 > T1,T2»YI 
IF(T1.  LE.  0.  ) GO  TO  15 
TT=T1  -YI 

I Z=3 

T=TTT  <T) 

304  .TT=TTi-YI 
T=TTT (TT) 

IFCTT.  GT.  T2>  GO  TO  302 
CALL  RB(T> 

GO  TO  27 
9 STOP 
END 


o o o 
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CORRESPONDING  STATES  CORRELATION  OF  THE  DIRECT 
CORRELATION  FUNCTION  INTEGRALS  OF  PURE  LIQUIDS 
MAIN  PROGRAM 
IMPLICIT  REAL*8(A-H, O-Z) 

DIMENSION  P ( 4 > , STEP ( 4 ) , AB<4>,  BB<4).  NAME (20) 

DIMENSION  TEMP ( 22 ) , TR(22), DR(8>, RRHO(S), DEXP(8), 

1DCAL ( S ) » DEV (8).  ERR (8) 

DIMENSION  VI (8, 22), Cl (8,  22), PI (8,  22) , PI  I (8), 

1PRES ( 8 ) , V<8),VII<8),CII(8> 

COMMON/DATA/TEMP,  VI,  Cl,  DEXP,  DCAL,  VC 
COMMON/COND/NP,  NU,  NN,  TERM 
COMMON/PR INT/DDEX (8,  22),  DDCA(8,  22) 

COMMON/PARAM/RHOS,  DCFS,  TS 
READ (5,1)  NAME 
1 FORMAT (20A4) 

READ (5,  3)  NJ 

3 FORMAT (12) 

WRITE (6,4)  NU 

4 FORMAT (IX, 'NO.  OF  ISOTHERMS=  ',12,//) 

READ(5,  20)  NP,  10 

20  FORMAT (212) 

READ( 5, 7)  NPASS, (P(I), 1=1, NP ) , (STEP( I), 1=1, NP ) 

READ (5,  8)  (AB(  I >,  1 = 1,  NP  > , (BB(  I ),  1 = 1,  NP) 

7 FORMAT( 12,  6F10.  0) 

8 FORMAT (3X,  6F10.  0) 

WRITE (6, 101 ) 

101  FORMAT (IX,  'PATTERN  SEARCH  ESTIMATION  OF  PARAMETERS  FROM 
1DCFI  DATA '/IX, 'CONTROL  PARAMETERS  WERE  SET  AS  FOLLOWS-'/) 

WR ITE (6,  102)  NP,  10 

102  FORMAT< IX, 5H  NP=,  I2/1X,  5H  I0=,  12) 

WRITE (6, 106)  NPASS 

106  FORMAT ( IX, 15HNUMBER  OF  PASS=, 12) 

WRITE (6,  107) 

107  FORMAT( IX, 18HINITIAL  PARAMETERS) 

WRITE(6,  108)  (I, P(  I),  1 = 1, NP) 

108  FORMAT( IX,  12, F10.  3) 

WR ITE( 6,  109) 

109  FORMAT( IX,  12HSIZE  OF  STEP) 

WRITE(6, 110)  (I, STEP(I), 1=1, NP) 

110  FORMAT( IX,  12, F10.  3) 

WRITE ( 6,  111 ) 

111  FORMAT (IX,  30HL0WER  LIMITS  OF  THE  PARAMETERS) 

WR ITE (6, 112)  ( I, AB < I ) , 1=1, NP ) 

112  FORMAT( IX,  12, F10.  3) 

WRITE (6, 113) 

113  FORMAT ( IX, 30HUPPER  LIMITS  OF  THE  PARAMETERS) 

WRITE(6, 1 14)  (I, BB(I>, 1=1, NP ) 

114  FORMAT(  IX,  12,  F10.  3) 

READ (5,  120)  TC,  VC 

120  F0RMAT(2F10.  0) 

DO  200  11=1, NU 
READ (5, 300)  TEMP ( I I ) 

READ(  5,  300)  VII,CII,PII 
300  FORMAT (8F 10.  0) 

DO  400  KK=1,  8 
VI(KK, I I ) =VI I ( KK ) 

C I ( KK,  1 1 ) =C 1 1 ( KK ) 

PI(KK,  1 1 ) =P 1 1 ( KK ) 

400  CONTINUE 
200  CONTINUE 
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CALL  PATERNCNP, P, STEP. NPASS, 2, AE, BB  ) 

CALL  DCF I (P. SUM) 

WRITE (6, 202)  NAME 

202  FORMAT ( /// /20H  THE  SUBSTANCE  IS  : 20A4//) 

WR ITE  < 6,  130  > TC,  VC,  TS,  RHOS,  DCFS 
130  FORMAT  ( //  ' TC=',F8.3,  ' K'/'  VC=', 

1F8.  3,  ' CC/GMOL'/'  T*=', 

1F8.  3,  ' 'A'  / ' D*=',F8.  4,  ' GMOL/L  ' / ' C*=',F9.  4//) 

T0T=0. ODO 

DO  116  J=l, NJ 

TR(U)=TEMP(J)/TS 

WR I TE  < 6,  1 1 7 ) TEMP  < J ) , TR  < J ) 

117  FORMAT (//IX,  'FOR  THE  ISOTHERM  AT',  F7.  2,  ' DEG.  K,  ', 

1 'THE  REDUCED  TEMPERATURE  IS',F7.  4//) 

WRITE<6, 115) 

115  FORMAT (2X, 3HN0.  , 3X, 5HP, BAR, 3X,  8HV, CC/MOL,  5X, 4HD/D*, 4X, 
15HC,  EXP,  5X,  5HC,  CAL,  7X,  4HDEV.  , 5X,  5H7.  ERR/) 

DO  113  1 = 1,8 
PRES ( I ) =P I ( I , J) 

V(I)=VI(I,  J) 

DEXP ( I ) =DDEX ( I , J) 

DCAL( I )=DDCA< I,  J) 

IF ( DABS ( DEXP  < I ) > . LT.  1. OD-5 ) GO  TO  118 
DEV ( I ) =DEXP ( I ) -DCAL ( I ) 

ERR ( I ) =DEV ( I ) /DEXP  < I ) *100.  DO 
TOT =TOT +DABS ( ERR ( I ) ) 

DR ( I >=VC/V( I > 

RRHO ( I ) =1 . 0D03/VI ( I,  J)/RHOS 

WRITE (6, 119)  I, PRES ( I ) , V( I ) , RRHO ( I ) , DEXP ( I ) , 

1DCALCI),  DEV ( I ) , ERR ( I ) 

119  FORMAT (IX,  13,  F10.  2,  7F10.  4 ) 

118  CONTINUE 

116  CONTINUE 
AVG=TOT /TERM 
WRITE<6,  201)  AVG 

201  FORMAT ( / /30H  AVERAGE  PERCENTAGE  ERROR  IS  =, F10.  4//) 
STOP 
END 

SUBROUTINE  PATERN ( NP, P, STEP, NPASS,  10, AB, BB) 

IMPLICIT  REAL*8 C A— H, Q-Z ) , INTEGER(I-N) 

DIMFNSION  P ( NP ) , STEP ( NP ) 

DIMENSION  Bl(4>, B2(4), T(4>, S ( 4 ) , AB ( 4 ) , 33(4) 

I OF  =6 

NRD=NPASS 

L=1 

ICK=2 

ITTER=0 

DO  5 1 = 1,  NP 

B 1 ( I ) =P ( I ) 

B2( I )=P( I ) 

T ( I > =P  < I ) 

5 S(  I >=STEP(  I )*10. 

CALL  BOUNDSCP, IOUT, AB, BB, NP) 

IF ( IOUT.  LE.  0)  GO  TO  10 
IF( 10.  LE. 0)  GO  TO  6 
WR ITE( I OF, 1005) 

WRITEdOF,  1000)  (J,  P(U>,  U=l,  NP  ) 

6 GO  TO  999 

10  CALL  DCFI (P, Cl ) 

IF(  10.  LE.  1 ) GO  TO  11 
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WRITE! IOF, 1001 > ITTER,C1 
WRITE (I OF i 1000)  !J,  P(J),  J=l,  NP) 

11  DO  9?  INRD=1 , NRD 
DO  12  1=1, NP 

12  S ! I ) =S ! I ) / 10. 

IF(  10.  LE.  1 ) GO  TO  20 
WR ITE( I OF, 1003) 

WRITE! IOF, 1000)  ( J, S ( J ) , U=1 , NP ) 

20  IFAIL=0 

DO  30  1=1, NP 
IC=0 

21  P! I )=T! I )+S!  I > 

IC=IC+1 

CALL  BOUNDS ( P, IOUT, AB, BB, NP) 

IF ( IOUT.  GT.  0 ) GO  TO  23 
CALL  DCFI (P, C2 ) 

L=L+1 

IF ( 10.  LT.  3)  GO  TO  22 
WRITE (I OF,  1002)  L,  C2 
WRITE! I OF, 1000)  ( J, P ( J) , J=l, NP ) 

22  IF ! C 1 -C2 ) 23,23,25 

23  IF( IC.  GE.  2)  GO  TO  24 
S(I)=-S(I) 

GO  TO  21 

24  IFAIL= IFAIL+1 
P! I >=T ! I ) 

GO  TO  30 

25  T ( I ) =P ( I ) 

C 1=C2 

30  CONTINUE 

IF( IFAIL.  LT.  NP)  GO  TO  35 
IF ( ICK.  EQ.  2)  GO  TO  90 
IF ( ICR.  EQ.  1)  GO  TO  35 
CALL  DCFI (T, C2) 

L=L+ 1 

IF ( 10.  LT.  3)  GO  TO  31 
WRITE! IOF,  1002)  L,  C2 
WRITE! IOF, 1000)  ! J, T! J) , J=l, NP > 

31  IF ! C 1-C2 ) 32,34,34 

32  ICK=1 

DO  33  1=1, NP 
B 1 ! I >=B2( I ) 

P ! I )=E2( I ) 

33  T( I )=B2! I ) 

GO  TO  20 

34  C1=C2 

35  IB  1=0 

DO  39  1=1, NP 
B2( I )=T! I ) 

IF! DABS! B1 ! I >-B2! I ) ) . LT.  < DABS ! S ! I ) )*. 01 ) > IB1  = IB1  + 1 

39  CONTINUE 

IF!  IB1.  EQ.  NP)  GO  TO  90 
ICK=0 

ITTER=ITTER+1 

IF!  10.  LT.  2)  GO  TO  40 

WRITE! IOF, 1001 ) ITTER,C1 

WRITE! IOF, 1000)  !U, T! J) , J=l, NP > 

40  SJ=1. 

DO  45  11=1, 11 
DO  42  1=1, NP 
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T ( I > =B2 ( I > +SJ* < B2 ( I ) -B 1 ( I ) ) 

42  P ( I ) =T( I > 

SJ=SJ-.  1 

CALL  BOUNDS ( T , I OUT , AB , 33 , NP ) 

IF ( I OUT.  LT.  1 > GO  TO  46 
IF < II.  EQ.  11  ) ICK=1 

45  CONTINUE 

46  DO  47  1=1, NP 

47  BKI  >=B2(  I > 

GO  TO  20 

90  DO  91  1=1, NP 

91  T ( I ) =B2 ( I ) 

99  CONTINUE 

DO  100  1=1, NP 
100  P< I >=7  < I ) 

COST=C 1 

IF ( 10.  LE.  0)  GO  TO  999 

WRITE( IOF,  1004)  L,  Cl 

WRITE( IOF, 1000)  (J,  P<  J),  J=l,  NP) 

999  CONTINUE 
RETURN 

1000  FORMAT( IX,  5( 17, E13.  6)/> 

1001  F0RMAT(//1X,  13HITERATI0N  NO.  , 1 5/5X,  5HC0ST=, El  5.  6, 20X, 

1 10HPARAMETERS ) 

1002  FORMAT ( 10X, 3HN0.  , 14, 8X, 5HC0ST=, El 5.  6) 

1003  FORMAT (/IX, 28HSTEP  SIZE  FOR  EACH  PARAMETER) 

1004  FORMAT (IX,  13HANSWERS  AFTER,  13,  2X,  22HFUNCTI0NAL  EVALUATIONS 
1//5X,  5HC0ST=, El 5.  6, 20X,  18H0PTIMAL  PARAMETERS) 

1005  FORMAT (IX, 35HINITIAL  PARAMETERS  OUT  OF  BOUNDS  ) 

END 

SUBROUTINE  BOUNDS(P, IOUT, AB, BB, NP ) 

REAL*8  AB ( 4 ) , BB ( 4 ) , P ( 4 ) 

I0UT=0 
DO  1 1=1, NP 

1 IF(P(  I).  LE.  AB(  I ) .OR.  P(  I).  GT.  BB(I>  ) I0UT=1 
RETURN 
END 

SUBROUTINE  DCFI ( P, SUM ) 

IMPLICIT  REAL*8(A-H, 0-Z) 

DIMENSION  VI (8, 22) , Cl (8, 22) , RH0(8) 

DIMENSION  P ( 4 ) , TEMP(22), DEXP(8), DCAL(8) 

COMMON/DATA/TEMP, VI, Cl, DEXP, DC AL, VC 
COMMON/COND/NP, NJ, NN, TERM 
COMMON/PR INT/DDEX( 8, 22), DDCA(8, 22) 

COMMON/PARAM/RHOS,  DCFS,  TS 
RHOS=P ( 1 ) 

DCFS=P ( 2 ) 

TS=P ( 3 ) 

SUM=0.  ODO 
TERM=0.  ODO 
DO  6 J=l,  NJ 
TT=TS/TEMP ( J) 

NN=0 

DO  10  K»l,  8 

IF(VI  (K,  J).  LT.  1.  ODO)  GO  TO  10 
NN=NN+1 

RH0(NN)  = 1.  0D03/VI ( K,  J) 

DEXP ( NN > =C I ( K,  J) 

10  CONTINUE 

DO  13  1=1, NN 
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DN=RHQ  ( I ) /RHQS 
TERM-TERM+1 

A0=0.  98642D01-0.  10191D02*TT-0.  1 5356D0 1 *TT*TT 
Al=-0. 28465D02-*-0.  30864D02*TT+0.  60294D01*TT*TT 
A2=0.  27542D02-0.  32898D02*TT-0.  87130D01«TT*TT 
A3=-0. 82606D01+0. 12737D02*TT+0. 40170D01*TT*TT 
DCAL( I ) = ( AO+DN*  < A1 +DN* ( A2+DN*A3 ) ) >*DCFS 
SUM=SUM+ ( DCAL ( I ) -DEXP ( I ) > **2 
DDEX ( I / J ) =DEXP ( I ) 

DDCA< I, J ) =DCAL ( I ) 

13  CONTINUE 
6 CONTINUE 

SUM=DSQRT< SUM/ < TERM- 1.  ODO) ) 

RETURN 

END 
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C CORRESPONDING  STATES  CORRELATION  OF  THE  VOLUMETRIC 

C PROPERTIES  OF  COMPRESSED  PURE  LIQUIDS. 

C MAIN  PROGRAM 

IMPLICIT  REAL*8 ( A-H»  Q-Z ) 

DIMENSION  P ( 4 ) ( STEP ( 4 ) , AB < 4 > , BB<4),  NAME (20) 

DIMENSION  TEMP ( 34 ) , TR(34>,  DR<8),  RRH0(8),  PEXP(8), 

1PCAL ( 8 ) , DEV ( 8 ) , ERR (8) 

DIMENSION  VI <8, 34), Cl (8,  34) , PI (8,  34) , PI  I (8) , 

1PRES ( 8 ) » V(8)<VII(8)i  C 1 1 ( 8 ) 

COMMON /DATA /TEMP , VI,  PI,  PEXP,  PCAL,  VC 
COMMQN/COND/NP,  NJ,  NN,  TERM 
COMMQN/PR INT /PPEX  < 8,  34),  PPCA ( 8,  34) 

COMMON/PARAM/RHOS,  DCFS,  TS 
READ < 5,1)  NAME 
1 FORMAT ( 20A4 ) 

READ (5,  3)  NU 

3 FORMAT (12) 

WRITE ( 6,  4 ) NU 

4 FORMAT (IX,  'NO.  OF  ISOTHERMS=  ',12,//) 

READ(5, 20)  NP,  10 

20  FORMAT (212) 

READ (5, 7)  NPASS,  ( P ( I ) , 1 = 1, NP),  ( STEP ( I ) , 1 = 1, NP ) 

READ (5,  8)  (AB(  I ),  1 = 1,  NP),  (BB(  I ),  1 = 1,  NP) 

7 FORMAT( 12,  6F10.  0) 

8 FORMAT (3X,  6F10.  0) 

WR ITE ( 6, 101 ) 

101  FORMAT (IX#  'PATTERN  SEARCH  ESTIMATION  OF  PARAMETERS  FROM 

1 PVTX  DATA '/IX,  'CONTROL  PARAMETERS  WERE  SET  AS  FOLLOWS-'/) 
WRITE (6,  102)  NP,  10 

102  FORMAT ( IX,  5H  NP=,  I2/1X,  5H  I0=,  12) 

WRITE (6,  106)  NPASS 

106  FORMAT ( IX, 15HNUMBER  OF  PASS=, 12) 

WRITE (6.  107) 

107  FORMAT< IX, 18HINITIAL  PARAMETERS) 

WRITE(6,  108)  ( I, P( I ),  1 = 1, NP) 

108  FORMAT (IX,  I2,F10.  3) 

WR IT  E ( 6,  109) 

109  FORMAT ( IX, 12HSIZE  OF  STEP) 

WRITE(6, 110)  (I, STEP(I), 1=1, NP) 

110  FORMAT( IX,  12, F10.  3) 

WRIT  E(6,  111 ) 

111  FORMAT ( IX,  30HL0WER  LIMITS  OF  THE  PARAMETERS) 

WRITE(6,  112)  (I,  AB(I),  1 = 1,  NP) 

112  FORMAT( IX,  12, F10.  3) 

WRITE ( 6,  113) 

113  FORMAT ( IX,  30HUPPER  LIMITS  OF  THE  PARAMETERS) 

WRITE(6,  114)  (I,  BB(I),  1 = 1,  NP  )' 

114  FORMAT( IX,  12, F10.  3) 

READ (5,  120)  TC,  VC 

120  FORMAT < 2F 1 0.  0) 

DO  200  1 1 = 1, NJ 
READ (5, 300)  TEMP ( I I ) 

READ  ( 5,  300 ) VII.CII.PII 
300  F0RMAT(8F10.  0) 

DO  400  KK=1,  8 
VI(KK, II)=VII(KK> 

CI(KK, I I ) =C I I ( KK ) 

PI(KK,  II)=PII(KK) 

400  CONTINUE 
200  CONTINUE 
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CALL  PATERN!NP, P,  STEP, NPASS, 2, AB, BB) 

CALL  DCF  I (P, SUM) 

WRITE (6, 202)  NAME 

202  FORMAT !////2lH  THE  SUBSTANCE  IS  : 20A4//) 

WRITE (6,  130)  TC, VC, TS, RHCS, DCFS 
130  FORMAT  ( //  ' TC=',F8.3,  ' K'/'  VC=', 

1F8.  3,  ' CC/GMOL'/'  T*=', 

1F3.  3,  K ' / ' D*=  ' , F8.  4,  ' GMOL/L  ' / ' C*=',F9.4//) 

T0T=0.  ODO 
DO  116  J= 1 , N J 
TR ( J ) =TEMP ( J ) /TS 
WRITE(6,  117)  TEMP ( J ) , TR ( J ) 

117  FORMAT (//IX,  'FOR  THE  ISOTHERM  AT',  F7.  2,  ' DEG.  K, 

1 'THE  REDUCED  TEMPERATURE  IS',F7.  4) 

WR  ITE  ( 6,  134  > PI!1,U),  VI(1,U) 

134  FORMAT (//IX,  'THE  REF.  PRESSURE  IS',F9.  4,  ' 

1 BAR,','  THE  REF.  MOLAR  VOLUME  IS',F9.  4,  ' CC/GMOL'//) 

WRITE (6,  115) 

115  FORMAT ( 2X, 3HN0.  , 3X,  5HP, BAR, 3X, 8HV, CC/MOL,  5X, 4HD/D*, 4X, 
15HP,  EXP,  5X,  5HP,  CAL,  7X,  4HDEV.  , 5X,  5H7.  ERR/) 

DO  118  1=2,8 
PRES ! I > =P I ( I , U> 

V< I >=VI ( I,  J) 

PEXP ( I ) =PPEX ( I,  U) 

PCAL(I)=PPCA!I, J) 

IF ( DABS ( PEXP ( I ) ) . LT. 1. OD-5 ) GO  TO  118 
DEV ! I ) =PEXP ! I ) -PCAL ( I ) 

ERR ( I ) =DEV ( I > /PEXP ( I ) *100. DO 
TOT =TQT +DABS ( ERR ( I ) ) 

DR ( I ) =VC/V ( I ) 

RRHO! I ) = 1.  0D03/VI ( I , J)/RHOS 

WRITE! 6,  119)  I,  PRES! I ) , V! I ),  RRHO ! I ) , PEXP ( I ) , 
lPCAL(I), DEV ! I ) , ERR! I ) 

119  FORMAT!  IX,  13,  F10.  2,' 7F10.  4) 

US  CONTINUE 

116  CONTINUE 
AVG=TOT /TERM 
WRITE16, 201)  AVG 

201  FORMAT ! / /30H  AVERAGE  PERCENTAGE  ERROR  IS  =, F10. 4//) 
STOP 
END 

SUBROUTINE  PATERN ! NP , P, STEP, NPASS, 10, AB, BB) 

IMPLICIT  REALMS ( A-H, O-Z), INTEGER! I-N) 

DIMENSION  P ! NP ) , STEP ! NP ) 

DIMENSION  B1 (4),  B2(4>,  T(4),  S!4),  A3 ! 4 ) , BB ! 4 ) 

I OF =6 

NRD=NPASS 

L=1 

ICK=2 

ITTER-0 

DO  5 1 = 1,  NP 

B 1 < I ) =P ! I > 

B2! I )=P! I ) 

T ( I ) =P  < I ) 

5 S ! I ) =STEP  ( I ) * 1 0. 

CALL  BOUNDS!P, IOUT, AB, BB, NP ) 

IF!  IOUT.  LE.  0)  GO  TO  10 
IF!  10.  LE.  0)  GO  TO  6 
WRITE! IOF,  1005) 

WRITE! IOF,  1000) !J,  P! J),  J=l,  NP ) 
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6 GO  TO  999 

10  CALL  DCF I (P» Cl  > 

I F ( ID.  LE.  1 > GO  TO  11 
WRITEdOF,  1001  > ITTER; Cl 
WRITE(  IOF,  1000)  <J,  P(J>,  J=1,NP> 

11  DO  79  INRD=1 , NRD 
DO  12  1=1, NP 

12  S<  I )=S(  D/10. 

IF(  10.  LE.  1 ) GO  TO  20 
WR ITE( IOF,  1003) 

WRITEdOF,  1000)  (J, S( J), J=l, NP ) 

20  IFA IL*0 

DO  30  1=1, NP 
IC=0 

21  P( I >=T( I )+S( I ) 

IC=IC+1 

CALL  BOUNDS(P, IOUT, AB, BB, NP) 

IF(  IOUT.  GT.  0)  GO  TO  23 
CALL  DCF I (P,  C2 ) 

L=L+1 

IF ( 10.  LT.  3)  Ga  TO  22 
WRITEdOF,  1002)  L,  C2 
WRITEdOF,  1000)  <J, P< J), J=l, NP) 

22  IF ( C 1 -C2 ) 23,23,25 

23  IF ( IC.  GE.  2)  GO  TO  24 
S< I )=-S( I ) 

GO  TO  21 

24  IFAIL=IFAIL+1 
P<I )«T( I) 

GO  TO  30 

25  T ( I ) =P ( I ) 

C 1 =C2 

30  CONTINUE 

IFdFAIL.  LT.  NP)  GO  TO  35 
IF ( ICK.  EQ.  2)  GO  TO  90 
IF < ICK.  EQ.  1 > GO  TO  35 
CALL  DCFICT,  C2) 

L=L+1 

IF( ID.  LT.  3)  GO  TO  31 
WRITE( IOF,  1002)  L,  C2 
WRITEdOF,  1000)  ( J, T( J), U=l,  NP ) 

31  IF ( C 1 -C2 ) 32,34,34 

32  ICK=1 

DO  33  1=1, NP 
B 1 ( I ) =B2 ( I ) 

P( I )=B2( I ) 

33  T ( I ) =B2 ( I ) 

GO  TO  20 

34  C 1=C2 

35  IB  1=0 

DO  39  1=1, NP 
B2 ( I ) =T ( I ) 

IF ( DABS <B1(I)-B2(I) ).  LT.  ( DABS ( S ( I ) ) *.  01 ) ) IB1  = IB1  + 1 
39  CONTINUE 

IF ( IB1.  EQ.  NP)  GO  TO  90 
ICK=0 

ITTER-ITTER+1 
IF(  10.  LT.  2)  GO  TO  40 
WRITEdOF,  1001 ) ITTER, Cl 
WRITEdOF,  1000)  (J,  T(J),  J=l,  NP  ) 
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40  SJ=t. 

DO  45  11=1, 11 
DO  42  1=1, NP 

T < I ) =82 ( I > +SJ* ( B2 ( I ) -B 1 < I ) ) 

42  P ( I ) =T ( I > 

SJ=SJ-.  1 

CALL  HOUNDS CT, IOUT,  AB, BB, NP ) 

IF  ( IOUT.  LT.  1 ) GO  TO  46 
IF(  II.  EQ.  11  ) ICK=1 

45  CONTINUE 

46  DO  47  1=1, NP 

47  Sl( I > =B2 ( I ) 

GO  TO  20 

90  DO  91  1 = 1,  NP 

91  T ( I ) =B2 ( I ) 

99  CONTINUE 

DO  100  1 = 1,  NP 
100  P ( I ) =T ( I ) 

COST=C 1 

IFdO.  LE.  0)  GO  TO  999 

WRITEdOF,  1004)  L, Cl 

WRITE< I OF,  1000)  < J,  P< J) , J=l, NP  ) 

999  CONTINUE 
RETURN 

1000  FORMATC IX, 5( 17, E13. 6) /) 

1001  F0RMAT(//1X,  13HITERATI0N  NO.  , I 5/5X,  5HC0ST=, El  5.  6, 20X, 

1 1 OHPAR AMETERS ) 

1002  FORMAT ( 10X,  3HN0.  , 14,  SX,  5HC0ST=,  E15  6) 

1003  FORMAT (/IX,  2SHSTEP  SIZE  FOR  EACH  PARAMETER) 

1004  FORMATdX,  13HANSWERS  AFTER,  13,  2X,  22HFUNCTI0NAL  EVALUATIONS 
1//5X,  5HC0ST=, El 5.  6, 20X,  1SH0PTIMAL  PARAMETERS) 

1005  FORMATdX,  35HINITIAL  PARAMETERS  OUT  OF  BOUNDS  ) 

END 

SUBROUTINE  BOUNDS(P, IOUT, AB, BB, NP ) 

REAL*8  AB ( 4 ) , BB ( 4 ) , P ( 4 ) 

IOUT  =0 
DO  1 1=1, NP 

1 IF(P(I).  LE.  ABU)  .OR.  P ( I ) . GT.  BB  ( I ) ) I0UT=1 
RETURN 
END 

SUBROUTINE  DCFI(P,SUM> 

IMPLICIT  REAL*8<A-H, 0-Z) 

DIMENSION  VI (8, 34), PI (8, 34),  RHO (8) 

DIMENSION  P ( 4 ) , TEMP ( 34 ) , PEXP  <8 ) , PCAL(8) 

COMMON/DATA/TEMP,  VI,  PI,  PEXP,  PCAL,  VC 
COMMON/COND/NP, NJ,  NN,  TERM 
COMMON/PRINT /PPEX ( S,  34),  PPCA<8,  34) 

COMMON/PARAM/RHOS,  DCFS,  TS 
RHQS=P ( 1 ) 

DCFS=P ( 2 ) 

TS=P (3) 

RK=0.  0831 434D0 
SUM=0.  ODO 
TERM=0.  ODO 
DO  6 J= 1 , N J 
TT=TS/TEMP(J) 

NN=0 

DO  10  K=1 , 8 

IF(VI(K,  J).  LT.  1.  ODO)  GO  TO  10 
NN=NN >• 1 
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RH0CNN>=1. 0D03/VI (K, J) 

DR=RHO< 1 ) /RHOS 
PEXP ( NN ) =P I (K, J) 

PR=PEXP ( 1 > 

10  CONTINUE 

DO  13  1=2, NN 
DN=RHO ( I ) /RHOS 
TERM=TERM+1 

A0=0.  98642D01— 0.  10191D02*TT-0.  15356D01*TT*TT 
Al=-0. 28465D02+0. 30864D02*TT+0. 60294D0 1 *TT*TT 
A2=0.  27542D02-0.  32898D02*TT-0.  37130D01*TT*TT 
A3=-0.  82606DG1+0.  12737D02*TT+0.  4G170D01*TT*TT 
BO=RHOS* ( DN-DR ) *RK*TS/TT 
B 1=DN  >-DR 

B2=DN*DN+DN*DR+DR*DR 
B3=DN*DN*DN+DN*DN*DR+DN*DR*DR+DR*DR*DR 
C0=A0>-Al*Bl/2.  0D0+A2*B2/3.  0D0+A3*B3/4.  ODO 
PCAL( I )=PR+BO*< 1.  ODO-DCFS*CO> 

SUM=SUM+(PCAL< I )-PEXP ( I ) ) **2 
PPEX(I,J)=PEXP(I) 

PPCA( I,  J)=PCAL( I ) 

13  CONTINUE 
6 CONTINUE 

SUM=DSQRT < SUM/ < TERM- 1.  ODO) ) 

RETURN 

END 


o o o 
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CORRESPONDING  STATES  CORRELATION  OF  THE  VOLUMETRIC 
PROPERTIES  OF  COMPRESSED  LIQUID  MIXTURES. 

MAIN  PROGRAM 

IMPLICIT  REALMS ( A-H, O-Z ) 

DIMENSION  P ( 4 ) , STEP ( 4 ) , AB < 4 ) , BB ( 4 ) , NU ( 8 ) , NAME (20) 
DIMENSION  TC(3>,  VC (3),  TS(3), VS (3), CS(3> 

DIMENSION  TCM ( 8 ) » VCM(8>,  TSM(8) , RH0S(8), DCFS (8) 

DIMENSION  XI(8, 3), TEMP(24, 8), TR(24, 8), RRH0(8) 

1 , PEXP ( 8 ) , PCAL ( 8 ) , DEV ( 8 > , ERR ( 8 ) 

DIMENSION  VI  (8,  24,  8),  Cl  (8,  24,  8),  PI  (8,  24,  8) 

1,  VII('3),  CII(S),  P 1 1 ( 8 ) , V(8>,  PRES  ( 8 ) 

CQMMON/DATA/TEMP, VI,  PI,  PEXP, PCAL 
COMMON/COND/NP, NC,  NI, NJ, NN, TERM 
COMMON/PR INT/PPEX( 8,  24, 8), PPCA<8, 24, 3) 

COMMON /PAR AM /TCM,  VCM,  TSM, RHOS,  DCFS 
COMMON/CRIT/XI, TC, VC, TS, VS, CS 
READ (5, 13)  NAME 
13  FORMAT ( 20A4 ) 

READ ( 5, 101 ) NP, 10 

101  FORMAT (212) 

READ (5,  102)  NPASS,  (P ( I ) , 1 = 1, NP  >,  ( STEP ( I ) , 1 = 1, NP) 

READ (5, 103)  ( AB ( I ) , 1=1, NP), (BB(I), 1=1, NP) 

102  FORMAT (12, 6F10.  0) 

103  FORMAT ( 3X, 6F 10.  0) 

WRITE(6, 104) 

104  FORMAT (IX, 'PATTERN  SEARCH  ESTIMATION  OF  PARAMETERS  FROM 
1DCFI  DATA '/IX,  'CONTROL  PARAMETERS  WERE  SET  AS  FOLLOWS-'/) 

WRITE(6, 105)  NP,  10 

105  FORMAT( IX,  5H  NP=,  I2/1X,  5H  I0=,  12) 

WRITE (6,  106)  NPASS 

106  FORMAT ( IX,  15HNUMBER  OF  PASS=,  12) 

WRITE (6,  107) 

107  FORMAT( IX, 18HINITIAL  PARAMETERS) 

WRITE(6,  108)  ( I, P(  I ),  1 = 1, NP ) 

108  FORMAT( IX,  12, F10. 3) 

WRITE (6,  109) 

109  FORMAT( IX, 12HSIZE  OF  STEP) 

WRITE(6,  110)  < I , STEP ( I > , I = 1 , NP  > 

110  FORMAT( IX,  12, F10.  3) 

WR ITE ( 6,  111 ) 

111  FORMAT (IX, 30HL0WER  LIMITS  OF  THE  PARAMETERS) 

WR ITE( 6,  1 12 ) (I, AB(I),  1 = 1, NP > 

112  FORMAT( IX,  12,  F10. 3) 

WRITE(6,  113) 

113  FORMAT (IX,  30HUPPER  LIMITS  OF  THE  PARAMETERS) 

WRITE(6,  114)  (I, BB(I),  1 = 1, NP ) 

114  FORMAT( IX,  12, F10.  3) 

READ (5,  1 ) NC 

1 FORMATX 12 ) 

WRITE (6,  2)  NC 

2 FORMAT  (/IX,  'NO.  OF  COMPONENTS=  ',  I2/>- 
READ  (5,3)  (TC(I),  VC  < I ) , 1 = 1,  NC  ) 

READ (5, 3)  <TS( I), VS (I ), CS( I ),  1 = 1, NC ) 

3 FORMAT (8F10.  0) 

READ ( 5,  1 ) NI 
DO  100  11  = 1,  NI 

READ(  5,  3)  (XKII,  I),  1 = 1,  NC) 

READ( 5, 1 ) NU( II ) 

DO  200  JJ=1, NU ( 1 1 ) 

READ (5,3)  TEMP ( JJ, II) 
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READ (5,  3)  VII,CII,PII 
DO  300  KK=1 , 8 
VI  (KK, JJ, II)=VII(KK) 

Cl (KK, JJ, II )=CII (KK> 

PI  (KK,  JJ,  II  >=PII (KK) 

300  CONTINUE 

200  CONTINUE 
100  CONTINUE 

CALL  PATERN ( NP , P,  STEP,  NPASS,  2,  A3,  BB> 

CALL  DC FI (P, SUM) 

WR  IT  E<  6,.  301  > NAME 

301  FORMAT (///17H  THE  SYSTEM  IS  : 20A4///) 

TOT=G.  ODO 

DO  115  11=1, NI 

WRITE (6,  8)  CXICII,  I),  1 = 1, NC-1) 

8 FORMAT ( ///IX,  ' X1=',2F9.  5) 

WR I TE < 6, 9 ) TCM (II), VCM (II), TSM (II), RHOS (II), DCFS (II) 

9 FORMAT ( ' TCM=',F8.  3,  ' K ' / ' VCM=', 

1FS.  3,  ' CC/GMOL'/'  T*=',F9.3,  ' K' 

1/'  D*=',F8.  4,  ' GMOL/L '/  ' C*=',F9.  4//) 

DO  116  JJ=1,  NJ( II ) 

TR ( JJ,  1 1 ) =TEMP ( J J,  1 1 ) /TSM (II) 

WRITE (6,  10)  TEMP ( JJ,  II), TR(JJ,  II) 

10  FORMAT (//IX,  'FOR  THE  ISOTHERM  AT',  F7.  2,  ' DEG.  K,  ', 

1 'THE  REDUCED  TEMPERATURE  IS',F7. 4) 

WRITt(6,  14)  PI(1, JJ,  II), VI <1,  JJ,  II) 

14  FORMAT < // 1 X,  'THE  REFERENCE  PRESSURE  IS',  F8.  4,  ' 

1 BAR,  ',  'THE  REFERENCE  MOLAR  VOLUME  IS',F9.  4//) 

WRI7E(6,  11) 

11  FORMAT ( 2X,  3HN0.  , 3X,  5HP,  BAR,  3X,  8HV,  CC/MOL,  5X,  4HD/D*,  4X, 
15HP,  EXP,  5X,  5HP,  CAL,-7X,  4HDEV.  , 5X,  5H7.  ERR/) 

DO  117  KK=2, 8 

PRES ( KK ) =P I ( KK,  JJ,  II) 

V ( KK ) —VI ( KK,  JJ,  II) 

PEXP(KK)=PPEX(KK,  J J,  1 1 ) 

PCAL < KK ) =PPCA ( KK,  JJ,  1 1 ) 

IF(DABS(PEXP (KK) ).  LT.  1.  OD-5)  GO  TO  117 
DEV ( KK ) =PEXP ( KK ) -PC AL ( KK ) 

ERR (KK ) =DEV ( KK ) /PEXP ( KK ) *100.  DO 
TOT=TOT+DABS ( ERR ( KK ) ) 

RRHO ( KK ) = 1 . 0D03/VI (KK,  JJ,  II )/RHOS( II ) 

WRITE (6,  12)  KK, PRES(KK) , V(KK> , RRHO(KK) , PEXP (KK) , 
IPCAL(KK) , DEV(KK) , ERR(KK) 

12  FORMAT(  IX,  13,  F10.  2,  7F10.  4) 

117  CONTINUE 

116  CONTINUE 
115  CONTINUE 

AVG=TQT/TERM 
WRITE( 6,  201 > AVG 

201  FORMAT (//30H  AVERAGE  PERCENTAGE  ERROR  IS  =, F10. 4//) 
STOP 

END 

SUBROUTINE  PATERN ( NP , P, STEP, NPASS,  10, AB, BB ) 

IMPLICIT  REAL*8(A-H,  0-Z),  INTEGER ( I-N ) 

DIMENSION  P(NP), STEP(NP) 

DIMENSION  Bl(4), B2 ( 4 ) , T(4), S(4), AB ( 4 ) , BB ( 4 ) 

I OF  =6 
NRD=NPASS 
L=1 
ICK=2 
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ITTER=0 
DO  5 1 = 1,  NP 
S 1 ( I ) =P  ! I ) 

E2! I )=P( I) 

T ( I ) =P ( I > 

5 S< I >=STEP< I>*10. 

CALL  BOUNDS!P,  IOUT, AB, BB, NP  ) 

I F ( I OUT . LE.  O)  GO  TO  10 
IF ( 10.  LE.  0)  GO  TO  6 
WRITE! IOF,  1005) 

WRITE! I OF,  1000) ! J,  P ! U > , J=1 , NP ) 

6 GO  TO  999 

10  CALL  DCFI (P, Cl ) 

IF!  10.  LE.  1 ) GO  TO  11 
WRITE! IOF, 1001 ) ITTER,C1 
WRITE! IOF,  1000)  ( J, P ! J ) , J=1 , NP ) 

11  DO  99  INRD=1, NRD 
DO  12  1=1, NP 

12  S(I)=3(I)/10. 

IF! 10.  LE.  1 ) GO  TO  20 
WRITE! IOF, 1003) 

WRITE! IOF, 1000)  ( J, S ! J ) , J=1 , NP ) 

20  IFAIL=0 

DO  30  1=1, NP 
IC=0 

21  P ( I ) =T ! I ) +S ! I ) 

IC=IC+1 

CALL  BOUNDS!P, IOUT, AB, 3B, NP ) 

IF! IOUT.  GT.  0)  GO  TO  23 
CALL  DCFI ! P,  C2 ) 

L=L  + 1 

IF! 10.  LT.  3)  GO  TO  22 
WRITE! IOF,  1002)  L,  C2 
WRITE! IOF,  1000)  ( J, P ! J ) , J=1 , NP > 

22  IF ! C 1-C2 ) 23,23,25 

23  IF! IC. GE.  2)  GO  TO  24 
S ( I ) = -S  ( I ) 

GO  TO  21 

24  IFAIl =IFAIL+1 
P ( I ) =T ( I ) 

GO  TO  30 

25  T ( I ) =P ( I ) 

C 1=C2 

30  CONTINUE 

IF! IFAIL.  LT.  NP ) GO  TO  35 
IF! ICK.  EQ.  2)  GO  TO  90 
IF!  ICK.  EQ.  1 ) GO  TO  35 
CALL  DCFI !T, C2> 

L=L+-1 

IF!  10.  LT.  3)  GO  TO  31 
WRITE! IOF, 1002)  L, C2 
WRITE! IOF, 1000)  ( J, T( J ) , J=1 , NP ) 

31  IF ! C 1-C2 ) 32,34,34 

32  I CK=1 

DO  33  1=1, NP 
B 1 ! I )=B2! I ) 

P ( I )=B2< I > 

33  T ! I > =B2 ! I ) 

GO  TO  20 

34  C 1=C2 
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35  IB  1=0 

DO  3?  1 = 1,  NP 
B2( I )=T( I > 

IF  ( DABS  ( B 1 ( I >-B2< I ) ) . LT.  ( DABS  ( S ( I ) >*.  01  ) ) IB1  = IE1  + 1 

39  CONTINUE 

IF(  IB1.  EQ.  NP  > GO  TO  90 
ICK=0 

ITTER=ITTER+1 

IF(  10.  LT.  2)  GO  TO  40 

WRITE! I OF,  1001 > ITTER,C1 

WRITfc ! IOF,  1000)  ! J,  T! J),  J=l, NP > 

40  SJ=1. 

DO  45  11=1, 11 
DO  42  1=1, NP 

T( I )=B2( I >+SU*!B2( I >-Bl ! I ) ) 

42  P ( I ) =T ( I ) 

SJ=SJ-.  1 

CALL  BOUNDS ( T, IOUT, AB, BB, NP) 

IF!  IOUT.  LT.  1 ) GO  TO  46 
IF!  II.  EQ.  11 ) ICK=1 

45  CONTINUE 

46  DO  47  1=1, NP 

47  B 1 < I ) =B2 ! I ) 

GO  TO  20 

90  DO  91  1 = 1,  NP 

91  T( I >=B2< I ) 

99  CONTINUE 

DO  100  1=1, NP 
100  P < I ) =T ! I > 

COST=-C  1 

IF!  10.  LE.  0)  GO  TO  999 

WRITE! IOF, 1004)  L, Cl 

WRITE! IOF, 1000)  ! J, P ! U ) , J=1 , NP ) 

999  CONTINUE 
RETURN 

1000  FORMAT! IX, 5! 17, E13.  6)/) 

1001  F0RMAT!//1X,  13HITERATI0N  NO.  , I5/5X, 5HC0ST=, E15.  6, 20X, 

1 10HPARAMETERS) 

1002  FORMAT! 10X, 3HN0.  , 14, SX,  5HC0ST=, E15.  6) 

1003  FORMAT! /IX, 28HSTEP  SIZE  FOR  EACH  PARAMETER) 

1004  FORMAT  (IX,  13HANSWERS  AFTER,  13, 2X,  22HFUNCTI0NAL  EVALUATIONS 
1//5X,  5HC0ST=, El 5.  6, 20X,  18H0PTIMAL  PARAMETERS) 

1005  FORMAT ( IX, 35HINITIAL  PARAMETERS  OUT  OF  BOUNDS  ) 

END 

SUBROUTINE  BOUNDS!P,  IOUT, AB,  BB, NP ) 

REAL*8  AB<4), SB (4), P(4> 

I0UT=0 
DO  1 1 = 1,  NP 

1 IF!P!  I ).  LE.  AB!  I > .OR.  P ! I ) . GT.  BB  ( I ) ) I0UT=1 
RETURN 
END 

SUBROUTINE  DCFI(P.SUM) 

IMPLICIT  REAL*8 ! A-H, 0-Z ) 

DIMENSION  TC ( 3) , VC (3),  TS!3>,  VS!3>,  CS(3),  NU ( 8 ) 

DIMENSION  TCM ! 8 ) , VCM!8) , TSM!8) , VSM(8) , RH0S!8) , DCFS(8) 
DIMENSION  VI! 8,  24,  8) , PI ( 8, 24,  8 > , DI J<3, 3 ) , RHO!S) 

DIMENSION  P!4),  XI (8,  3), TEMP! 24, 8), PEXP!8), PCAL!3> 
COMMON/DATA/TEMP,  VI,  PI, PEXP, PCAL 
COMMON/COND/NP, NC, NI, NU, NN, TERM 
COMMON /PR  I NT /PPEX ! 8, 24, 8), PPCA!8, 24, 8) 
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COMMGN/PARAM/TCM,  VCM, TSM, RHOS, DCFS 
COMMON /CR IT/X 1/ TC, VC, TS, VS, CS 
DIJ(1, 2 ) =P ( 1 ) 

DI J(2, 1 )=DIJ( 1, 2) 

DIJ( 1,  1 )=0.  DO 
DIJ<2,  2 ) =0.  DO 
SUM=0.  DO 
TERM=0. DO 
DO  1 11=1, NI 
TCM (II) =0.  DO 
VCM (II) =0. DO 
TSM (II) =0. DO 
VSM (II) =0.  DO 
DCFS  (II) =0.  DO 
DO  2 1 = 1,  NC 

TCM ( 1 1 ) =TCM< 1 1 ) +X I ( 1 1 , I)*TC(I> 

TSM (II) =TSM ( 1 1 ) +X I ( 1 1 , I)*TS(I) 

2 DCFS (II) =DCFS (II) +X 1(11/  I)*CS(I) 

DO  3 1 = 1,  NC 

DO  3 J=l, NC 

VCM ( 1 1) =VCM ( 1 1 ) +X I < 1 1 , I ) *X I ( 1 1 , V ) * ( ( VC ( I) +VC < J ) ) /2.  DO ) 

3 VSM (II) =VSM (II) +X I ( 1 1 , I)*XI (II, J)*( ( (VS(I)+VS(J) )/2.  DO) 
1*(  1. DO-DI J( I,  J) ) ) 

RK=0.  0831434D0 
DO  6 JJ=1, NJ( II > 

TT=TSM (II) /TEMP ( JJ,  II) 

NN=Q 

DO  10  KK=1 , S 

IF(VI(KK,  JJ,  II).  LT.  1.  ODO)  GO  TO  10 
NN=NN+ 1 

RHO ( NN  > = 1 . 0D03/VI (KK, JJ,  II ) 

RHOS ( 1 1 ) = 1 . 0D03/VSM (II) 

DR=RHO( 1 )/RHOS( II ) 

PEXP(NN)=PI(KK,  JJ,  II) 

PR=PEXP( 1 ) 

10  CONTINUE 

DO  13  KK=2, NN 
DN=RHQ ( KK ) /RHOS (II) 

TERM=TERM+1 

A0=0. 98642D01-0.  10191D02*TT-0.  1 5356D01*TT*TT 
Al=-0.  28465D02+0.  3C864D02*TT+0. 60294D01*TT*TT 
A2=0. 27542D02-0.  32898D02*TT-0.  87130D01*TT*TT 
A3=— 0. 82606D01+0.  12737D02*TT+0.  40170D01*TT*TT 
BO=RHOS ( 1 1 ) * ( DN-DR ) *RK*TSM  < 1 1 ) /TT 
B1=DN>-DR 

32=DN*DN+DN*DR+DR*DR 

B3=DN*DN*DN+DN*DN*DR+DN*DR*DR+DR*DR*DR 
C0=A0i-Al*Bl/2.  GD0+A2*B2/3.  0D0+A3*B3/4.  ODO 
PC AL ( KK ) =PR+BO* ( 1 . ODO-DCFS ( 1 1 )*C0) 

SUM=SUM+ ( PCAL ( KK ) -PEXP ( KK ) ) **2 
PPEX ( KK, JJ, II) =PEXP ( KK ) 

PPC A ( KK, JJ, II) =PCAL ( KK ) 

13  CONTINUE 
6 CONTINUE 
1 CONTINUE 

SUM=DSQRT( SUM/ ( TERM- 1.  ODO) ) 

RETURN 

END 
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